Graduate  Aeronautical  Laboratories 
California  Institute  of  Technology 

Pasadena,  California  91125 


Chemical  Reactions  in  Turbulent  Mixing  Flows 


Paiil  E.  Dimotakis*  and  Anthony  Leonard** 


Air  Force  Office  of  Scientific  Research 
Grant  No.  F49620-94-1-0353 
15  June  1994  through  28  February  1998 

Final  Report 


Issued:  28  May  1998 
Corrected:  15  September  1998 


*  John  K.  Northrop  Professor  of  Aeronautics  and  Professor  of  Applied  Physics 
**  Professor,  Aeronautics 


jjjig  qoALfn  ^ 


Approved  rolcaf'fs 

distribution  unlimited,  ^ 


19981020  096 


REPORTbbCUMENTATION  PAGE 


AFRL-SR-BL-TR-98- 


Public  reporting  burden  for  this  collection  of  Information  Is  estimated  to  average  1  hour  per  response,  including  the  tir  r  y  X  3  3 

the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  su 

reducing  tills  burden  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  121  ’  C*  ^  i  t( 

Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503  _  - _ 


1 ,  AGENCY  USE  ONLY  (Leave  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

blank)  28May98  Final  Technical  Report  6/15/94-2/28/98 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS 

(U)  Chemical  Reactions  in  Turbulent  Mixing  Flows  PE  -  61102F 

PR  -  2308 
SA  -  BS 


3  and  maintaining 
suggestions  for 
d  to  the  Office  of 


6.  AUTHOR(S) 

Paul  E.  Dimotakis  and  Anthony  Leonard 


G-F49620-94-.1-0353 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Graduate  Aeronautical  Laboratories 
California  Institute  of  Technology 
Pasadena,  C A  91125 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research 
1 10  Duncan  Avenue,  B 1 15 
Bolling  AFB  DC  20332-8050 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Unrestricted 


12b.  DISTRIBUTION  CODE 


Approved  for  public  relaas©  J 

distribution 


13.  ABSTRACT  ('Max/mu/n  200  Words; 

This  program  focused  on  fundamental  investigations  of  mixing,  chemical-reaction,  and  combustion 
processes,  in  turbulent,  subsonic ,.  and  supersonic  free-shear  flows.  The  program  was  comprised  of  an 
experimental  effort;  an  analytical,  modeling,  and  computational  effort;  and  a  diagnostics, 
instrumentation,  and  da ta-acquis it ion- development  effort,  with  significant  progress  in  each.  With 
regard  to  gas-phase  shear- layer  mixing  and  combustion,  effects  of  inf low/ initial  conditions, 
compressibility,  and  Reynolds  number  were  experimentally  investigated  and,  to  a  large  extent, 
clarified.  New  measures  to  characterize  level  sets  in  turbulence  were  developed  and  successfully 
employed  to  characterize  experimental  data  of  liquid-phase  turbulent- jet  flows  as  well  as  three- 
dimensional  direct-numerical-simulation  data  of  Rayleigh-Taylor- instability  flows.  The 
computational  effort  has  added  to  our  understanding  of  the  (H2+NO) /F2  chemical  system  employed  in 
the  shear- layer -mixing  investigations  as  well  as  mixing  in  high-speed  flows,  along  with  further 
developments  in  Riemann-Invariant -Manifold  gasdynamic  simulation  techniques  and  their  application  to 
unsteady  detonation  phenomena.  On  the  diagnostic  front,  developments  in  digital  imaging  and  Image 
Correlation  Velocimetry  have  continued,  and  been  used  to  investigate  turbulent- jet  mixing,  the 
unsteady  flow  over  an  accelerating  airfoil,  to  mitigate  aliasing  problems  in  the  computer 
reconstruction  of  (2+1) -dimensional  isosurface  data,  and  in  other  applications. 


14.  SUBJECT  TERMS 

Turbulence,  chemical -reacting  flows,  combustion,  supersonic  flow, 
Detonations,  hypersonic  propulsion,  level  sets 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 


15.  NUMBER  OF  PAGES 

71 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


1.  Introduction 


The  documentation  below  describes  work  on,  “Chemical  Reactions  in  Turbulent 
Mixing  Flows,”  performed  under  sponsorship  of  AFOSR  Grant  No.  F49620-94-1- 
0353.  This  work  was  focused  on  fundamental  investigations  of  mixing,  chemical- 
reaction,  and  combustion  processes,  in  turbulent,  subsonic,  and  supersonic  free- 
shear  flows.  The  program  was  comprised  of 

•  an  experimental  effort,  in 

0  chemically-reacting  and  non-reacting,  subsonic  and  supersonic  turbu¬ 
lent  free-shear  layers; 

o  liquid-  and  gas-phase  turbulent  jets  discharging  in  an  otherwise  qui¬ 
escent  reservoir, 

o  preliminary  work  on  transverse  jets  in  a  cross-flow,  and 
0  a  study  of  the  flow  field  over  accelerating  airfoils; 

•  an  analytical,  modeling,  and  computational  effort,  focusing  on  the  behavior 
of  hyperbolic  systems  with  strong  fronts  (shocks/detonations); 

and 

•  a  diagnostics,  instrtunentation,  and  data-acquisition-development  effort. 


Parts  of  this  effort  were  cosponsored  by  AFOSR  URI  Grant  No.  F49620-93-1- 
0338,  on  “Interaction  of  Chemistry,  Turbulence,  and  Shock  Waves  in  Hypervelocity 
Flow,”  with  some  instrumentation  support  derived  from  the  imaging  efforts  sup¬ 
ported  under  AFOSR/DURIP  Grant  No.  F49620-95-1-0199.  Some  experiments 
were  also  xmdertaken  as  part  of  work  on,  “Whole-field  measurements  of  turbulent 
flow  for  the  study  of  aero-optical  effects,”  under  AFOSR  Grant  No.  F49620-94-1— 
0283. 

The  work  outlined  above  is  documented  in  the  report  that  follows,  with  ad¬ 
ditional  detailed  docTunentation  in  the  publications,  theses,  and  reports  that  were 
written  in  the  course  of  this  research  program.  These  are  listed  in  Sec.  10,  below. 


2 


2.  Shear-layer  mixing  and  combustion 

The  experimental  data  discussed  in  this  section  were  investigated  in  the  GAL- 
CIT  Supersonic  Shear-Layer  Facility.  This  unique  facility  allows  the  study  of  high- 
Reynolds-number,  chemically-reacting,  turbulent  shear-layer  flows,  by  employing 
the  (H2  -I-  N0)/F2  chemical  system  for  “flip  experiments”  (Koochesfahani  &  Di- 
motahis  1986),  in  addition  to  color-schlieren  and  Rayleigh-laser  imaging,  for  flow 
visualization.  This  effort  has  focused  on  several  aspects  of  turbulent  shear-layer  be¬ 
havior,  spanning  incompressible-  to  compressible-flow  regimes.  In  particular,  the  re¬ 
sults  of  studies  on  effects  of  inflow/initial  conditions,  compressibility,  and  Reynolds 
number,  will  be  described  in  what  follows. 

As  discussed  by  various  authors,  e.g.,  Dimotakis  (1986),  high-Reynolds-number 
mixing  processes  can  be  conceptually  broken  into  three  separate  stages:  large-scale 
entrainment,  intermediate- scale  stirring,  and  small-scale  diffusion/mixing,  with  each 
stage  important  to  the  overall  molecular  mixing,  and  subsequent  combustion,  of  two, 
initially-separated,  freestream  fluids. 

This  sequence  can  be  expressed  in  terms  of  the  quantities  accessible  to  current 
measurement  techniques,  i.e., 


where  Sra  is  the  equivalent  mixed-fluid  thickness,  and  6  the  outer-scale  thickness.  In 
these  experiments,  St  is  employed  as  a  shear-layer  thickness  measure,  defined  as  the 
1%  temperature-rise  thickness,  i.e.,  the  transverse  extent  of  the  shear  layer  which 
lies  within  1%  of  the  maximiun  temperature  rise  (Dimotakis  1991a).  As  such,  Sm/Sx 
represents  the  fraction  of  molecularly-mixed  fluid  in  the  layer,  or  equivalently,  the 
molecular-mixing  efficiency.  The  first  term  represents  the  first  stage  of  mixing,  while 
the  second  represents  the  process  at  all  of  the  remaining  spatial  scales,  down  to  the 
diffusion,  or  Batchelor,  scale.  Because  of  the  large  dynamic  range  of  spatial  scales 
present  in  a  high-Reynolds-number  flow,  Sj^/S  is  inaccessible  to  direct-measurement 
techniques,  e.g.,  laser-induced  fluorescence,  and  consequently,  reliable  estimates 
of  molecular  mixing  must  rely  on  techniques  such  as  the  chemically-reacting  flip 
experiments  employed  here. 


AT/AT, 
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2.1  Inflow/initial  conditions 

A  set  of  experiments  was  performed  to  explore  effects  of  inflow  conditions,  at 
incompressible-flow  conditions.  At  a  (local)  Reynolds  number  of, 

_  pAU6(x) 


Reg  = 


~  2  X  10% 


and  freestream-velocity  and  -density  ratios  of, 

r  =  —  0-4  and  s  =  />2//>i— 1, 


respectively,  the  data  show  a  dramatic  influence  of  perturbed  inflow  conditions  on 
scalar  mixing. 


Fig.  1  Normalized  temperature-rise  data.  Left:  untripped  boundary  layers.  Right: 
tripped  high-speed  boundary  layer.  Reactant  compositions  -  diamonds:  <f)  = 
8;  triangles:  <f>  =  1/8-,  asterisks:  <{)  =  1/8,  reduced  chemical-kinetic  rate. 

Two  sets  of  temperature-rise  data  from  such  experiments  are  presented,  in 
Fig.  1,  with  corresponding  color-schlieren  flow- visualization  data  in  Fig.  2.  The 
temperature-rise  data,  in  the  form  of  the  “flip  experiment”,  provide  a  resolution-free 
probe  of  the  molecular-scale  mixing,  in  the  fast-kinetic  regime. 

The  only  difference  between  these  is  the  addition  of  a  0.8  mm-diameter  trip  wire 
on  the  splitter  plate  (high-speed  side),  50  mm  upstream  of  separation.  The  data 
are  measured  far  downstream  of  separation  (Bradshaw  1966),  x/9i  ~  3300,  where 
X  is  the  streamwise  coordinate  and  di  is  the  high-speed  boundary-layer  momentum 
thickness,  and  at  a  large  value  of  the  pairing  parameter  (Ho  L  Huang  1982  and 
Karasso  &  Mungal  1996),  P  ci  47.  These  measures  indicate  that  the  turbulence  can 
be  regarded  as  fully-developed. 
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Fig.  2  Color  schlieren  visualizations  of  untripped  (top)  and  tripped  (bottom)  flows. 

A  qualitative  difference  in  visualized-flow  structure  is  evident. 

This  change  in  inflow  conditions  has  a  significant  effect  on  all  scales  of  the 
flow:  the  shear-layer  growth  rate,  St/x,  has  decreased  by  21%;  the  mixed-fluid 
fraction  has  increased  by  11%;  and  the  mixed-fluid  composition  ratio  has 

decreased  by  9%.  Additionally,  the  data  indicate  a  change  from  a  non-marching 
scalar  probability-density  fimction  (pdf),  to  a  marching  pdf,  when  the  boundary 
layer  was  tripped.  These  observations  suggest  a  shear-layer  behavior  that  depends 
on  not  only  local-flow  properties,  but  also  on  upstream  conditions,  reminiscent  of 
results  from  low-dimensionality,  non-linear  (chaotic)  systems. 

A  more  complete  description  of  this  work  can  be  found  in  Slessor  et  al.  (1998a). 


2.2  Compressibility 

Compressibility  effects  on  turbulent  shear  layers  are  often  presented  in  terms 
of  the  total  convective  Mach  number  (Papamoschou  1989), 


Me  = 


U1-U2 

(Xi  0,2 


(2) 
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where  Ui  and  Cj  are  the  freestream  velocities  (speeds)  and  speeds  of  sound,  re¬ 
spectively.  In  the  following,  we  evaluate  the  utility  of  this  parameter,  as  well  aa 
investigate  the  effect  of  compressibihty  on  both  molecular  mixing  and  flow  struc¬ 
ture. 


2.2.1  Compressibility  and  growth 


Me 


Fig.  3  Normalized  compressible  shear-layer  growth  rates  as  parameterized  by  the 
total  convective  Mach  number.  Me  (Eq.  2). 


It  is  a  well-known  experimental  fact  that  compressibility  acts  to  reduce  to 
outer-scale  growth  rate,  6(x)/x.  In  particular,  compressibility  is  usually  assumed 
to  act  independently  of  freestream-velocity  and  -density  ratios,  in  the  form. 


6/x 

Sq/x 


_  U2  _  P2\ 

Pi) 


? 


(3) 
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with  estimates  for  the  incompressible-flow  growth  rate,  8q/x,  used  to  normaUze 
the  (measured)  compressible-flow  growth  rate,  6/x,  provided  by  the  spatial-growth 
model  of  Dimotahis  (1986).  Existing  compressible-flow  data  are  plotted  in  this 
fashion  in  Fig.  3. 

1.2 


1.0 


0.8 

\  0.6 
<0 

0.4 


0.2 


0.0 

0  12  3  4 

Fig.  4  Normalized  (as  in  Fig.  3)  compressible  shear-layer  growth  rates  parameter¬ 
ized  by  the  proposed  compressibility  parameter  He  (Eq.  4).  Solid  curve  is 
Eq.5. 

Abnormally-low  growth  rates  are  evident  in  Fig.  3,  and  suggest  either  an  inac¬ 
curate  representation  of  either  the  incompressible  shear-layer  growth  rate  or  the  flow 
compressibility.  Interestingly,  as  observed  by  several  authors,  e.ff.,  Hall  et  al.  (1993) 
and  Lu  &  Lele  (1994),  these  scaling  violations  are  observed  in  flows  with  extreme 
values  of  the  freestream  density  ratio,  s  =  P2/P1.  Existing  incompressible  growth- 
rate  data  (Brown  &  Roshko  1974)  are  well-represented  by  the  Dimotakis  (1986) 
model,  indicating  that  Mq  provides  an  inadequate  measure  of  compressibihty. 


AT/AT, 
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Reactant  compositions  -  diamonds:  (/>  =  8]  triangles:  =  1/8;  asterisks: 

^  =  1/8,  reduced  chemical-kinetic  rate. 

A  new  (maximum)  compressibility  scale  has  been  proposed, 


He 


max 

i 


Vii  - 1 

Oi 


X  AU, 


(4) 


derived  from  considerations  of  compressibility  as  a  kinetic-to-thermal-energy  con¬ 
version.  The  experimentally-measured  growth  rates,  normalized  as  in  Fig.  3,  are 
plotted  against  He  in  Fig.  4.  These  data  observe  a  near-functional  dependence  on 
this  parameter,  to  a  better  extent  than  on  Me,  with  a  solid  curve  through  the  data, 

~  (l -f  ,  Q!~4,  ;d~0.5.  (5) 

oo 

The  systematic  deviations  observed  in  flows  with  extreme  density/speed-of-sound 
ratios  are  absent  when  scaled  in  this  fashion. 


A  more  complete  discussion  of  these  results  is  available  in  Slessor  et  al.  (1998b). 


2.2.2  Compressibility  and  mixing 

A  set  of  experiments  was  performed  to  explore  effects  of  compressibility  on 
molecular  mixing.  The  partition  in  Eq.  1  indicates  that  compressibility  should  result 
in  a  net  decrease  in  the  amount  of  molecularly-mixed  fluid,  at  least  at  moderate- 
to-high  compressibility  conditions.  This  expectation  stems  from  two  experimental 
facts  (Dimotakis  1991a):  the  large-scale  growth  rate  is  approximately  1/2  of  its 
incompressible  value  by  Me  0.5,  and  the  mixed-fluid  fraction  is  approximately 
1/2  at  incompressible-flow  conditions. 
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At  nominally-constant  values  of  the  (local)  Reynolds  number, 

Re,  s  ~  9X10=, 

y- 

and  freestream-velocity  and  -density  ratios, 

r  =  1/2/171^0.4  and  s  =  /?2//0i^l, 
respectively,  flip  experiments  were  performed  at  two  compressibility  levels, 

Me  ~  0.25  and  Me  0.47 . 

Temperature-rise  data  from  these  two  chemically-reacting  flip  experiments  are  pre¬ 
sented  in  Fig.  5. 


Fig.  6  Color-schlieren  visualizations  of  low-compressibility  (Me  ~  0.25,  top)  and 
high-compressibility  (Me  ~  0.47,  bottom)  flows. 

Several  differences,  attributable  solely  to  the  increase  in  flow  compressibility, 
are  noted.  In  particular,  the  shear-layer  is  seen  to  be  thinner,  i.e.,  6/x  has  been 
reduced,  at  the  higher  compressibility,  in  quantitative  agreement  with  previous 
investigations  (Fig.  3).  For  the  H2-rich  {<f>  =  1/8)  parts  of  the  flip  experiment,  the 
maximum  normalized  temperature  rises  are  comparable  in  both  low-  and  moderate- 
compressibility  flows.  However,  the  maximum  attained  in  the  F2-rich  (^  =  8)  part 
is  significantly  (25%)  higher  than  in  low-compressibility  flow,  reaching  nearly  70% 
of  the  adiabatic  flame  temperature  for  this  mixture. 
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These  data,  despite  the  significant  increase  in  normalized  temperature-rise  at 
<f>  =  8,  show  only  a  modest  increase  in  mixed-fluid  fraction,  associated  with 

the  increase  in  Me.  This  reflects  a  change  in  the  shape  of  the  temperature-rise 
profile,  and  consequently,  the  mixed-fluid  structure.  Changes  are  also  observed  in 
the  instantaneous  spanwise-averaged  flow  visualizations  in  Fig.  6,  where  qualitative 
differences  in  visualized  structure  are  evident.  A  modified  flow  structure  is  also  in¬ 
dicated  by  the  change  in  the  mixed-fluid  composition  ratio,  which  is  seen  to  increase 
for  moderate-compressibility  flow,  favoring  the  high-speed  fluid  to  a  greater  extent 
than  in  the  corresponding  low-compressibility  flow. 

A  more  complete  description  of  these  experiments  is  available  in  Slessor  (1998). 


2.2.3  Compressiblity  and  flow  structure 

In  order  to  visualize  the  instantaneous,  spanwise-resolved  evolution  of  the  flow 
structure  with  increasing  compressibility,  laser-Rayleigh  scattering  images  were 
recorded  in  the  midspan  plane  of  three  different  shear-layer  flows.  These  data 
complement  both  the  spanwise-integrated  color-schlieren  data  and  flip-experiment 
temperature-rise  data  described  above.  The  images  were  recorded  with  a  pulsed 
(20  ns  FWHM  duration)  Nd:YAG  laser-illumination  source  formed  into  a  light  sheet, 
and  a  cryogenically-cooled  CCD  camera.  The  imaging  contrast  is  provided  by  the 
larger  index  of  refraction  (scattering  cross-section)  of  C2H4,  which  is  used  as  the 
low-speed-freestream  fluid,  in  conjunction  with  the  small  index  of  refraction  of  either 
N2,  or  He,  used  as  the  high-speed  freestream  fluid. 

Three  sample  images.  Figs.  7-9,  recorded  for  flows  at  compressibility  levels  of 
Me  0.15,  Me  ~  0.54,  and  Me  0.96,  respectively,  exhibit  a  qualitative  change 
with  increasing  compressibility.  In  particular,  the  low-compressibility  flow  is  charac¬ 
terized  by  typical  large-scale,  organized  “coherent  structure”,  which  is  less  evident 
at  mo  derate- compressibility  conditions,  and  even  less  so  at  high-compressibility  con¬ 
ditions.  These  findings  are  in  qualitative  accord  with  previous  investigations,  e.g., 
Clemens  k  Mungal  (1995).  In  the  highest-compressibility  flow,  a  traveling  shock- 
wave  system  has  been  captmed  in  the  lower,  low-speed  freestream,  and  is  similar 
to  those  observed  in  spanwise-integrated  (schlieren)  images,  at  similar  flow  condi¬ 
tions  (Hall  et  al  1993).  This  system  is  indicative  of  a  well-defined  flow  structure 
convecting  downstream  supersonically  with  respect  to  the  low-speed  freestream. 


Fig.  7  Rayleigh-scattering  visualization  of  a  subsonic  shear-layer  flow  at  ^  0.15. 
Flow  is  left-to-right,  with  the  high-speed  stream  on  the  top. 


Fig.  8  Rayleigh-scattering  visualization  of  a  moderate-compressibility  shear-layer 
flow  (Ml  ~  1.5 N2,  Me  ~  0.54). 
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Fig.  9  Rayleigh-scattering  visualization  of  a  high-compressibility  shear-layer  flow 
(Ml  ~  1.5  He,  Me  0.96). 

Even  though  a  reduction  in  organization  of  the  scalar  fleld  is  observed  with 
increasing  compressibility,  we  note  the  presence  of  distinct  scalar-fleld  interfaces 
in  all  images.  The  presence  of  such  interfacial  features  is  the  signature  of  large- 
scale  entrainment  and  mixing,  in  which  inducted  fluid  is  nearly  homogenized  by 
the  turbulent  stirring  and  diffusion  processes  (Dimotakis  1991a),  within  boundaries 
(separatrices)  across  which  transport  is  impeded.  As  evident  from  the  image  data 
presented  above,  these  interfaces  can  be  clearly  seen  in  all  three  compressibility 
cases  investigated. 

Further  experimental  details  and  a  discussion  of  these  investigations  in  the 
context  of  optical-beam  propagation  and  aero-optical  phenomena  can  be  found  in 
Dimotakis  et  al.  (1998b),  which  documents  parts  of  our  work  in  aero-optics,  “Whole- 
field  measurements  of  turbulent  flow  for  the  study  of  aero-optical  effects”.t 


t  Cofunded  by  AFOSR  Grant  No.  F49620-94-1-0283. 
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2.3  Reynolds-number  effects 

A  part  of  this  effort  has  been  to  document  the  effects  of  viscous  diffusion  on 
molecular  mixing,  at  high  Reynolds  numbers,  i.e..  Res  >  10^  (Dimotakis  1993). 
Although  these  effects  are  expected  to  be  small  over  the  range  of  Reynolds  numbers 
attainable  in  the  laboratory,  the  operating  range  of  practical  engineering  devices 
spans  many  orders  of  magnitude  (Dimotakis  1991a)  and  even  a  slow,  e.g.,  loga¬ 
rithmic,  dependence  on  Reynolds  number  can  amount  to  large  differences  between 
laboratory-scale  and  engineering-scale  devices. 

Open  symbols  -  natural  flows 


logio(Rex) 

Fig.  10  Shear-layer  6t  thickness  (cf.  Sec.  2)  data,  for  r  =  U2/U1  ~  0.4  and  s  = 
P2I pi  —  1,  plotted  as  a  function  of  (the  decimal  logarithm  of)  Reynolds 
number. 


Shear-layer  growth-rate  data,  recorded  for  several  different  chemically-reacting 
flows,  are  shown  in  Fig.  10,  plotted  vs. 


R.e^  — 


pUix 


(6) 


where  x  is  the  local  measuring  station,  to  uncouple  changes  in  the  local  shear- 
layer  thickness,  8{x),  from  other  changes  in  the  flow.  The  difficulty  of  ascertaining 
isolated  Reynolds  numbers  effects  is  evident  in  this  plot.  For  example,  the  effect 
of  tripping  the  high-speed  boundary  layer  (open/filled  symbols  joined  by  vertical 
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lines),  is  seen  to  have  a  more  significant  effect  on  the  flow  than  a  change  of  nearly 
two  orders  of  magnitude  in  Rcx- 

With  this  caveat,  the  shear-layer  mixed-fluid  fraction  data  are  presented  in 
Fig.  11,  along  with  predictions  by  models  of  Broadwell,  Breidenthal,  and  Mungal 
(BBM:  Broadwell  &  Mungal  1991);  and  of  Dimotakis  (1987).  The  inflow-condition 
experiments  described  above  axe  joined  by  a  line  (tripped  flow:  filled  symbol)  to 
indicate  the  magnitude  of  these  effects.  Results  from  Island  (1997)  are  also  included, 
along  with  an  estimate  for  the  renormalization  required  to  reconcile  their  choice  of 
a  5%-scalar-threshold  thickness  employed  in  normalizing  cold-chemistry  data,  with 
the  1%-temperature-rise  thickness,  St,  employed  in  normalizing  chemically-reacting 
data  in  the  work  undertaken  eis  part  of  this  Grant. 

It  is  perhaps  difficult  to  compare  the  validity  of  the  two  models,  on  the  basis 
of  these  data  alone.  However,  both  models,  as  well  as  the  data,  agree  with  an 
overall  downward  trend  of  the  mixed-fluid  fraction,  S^ISt,  with  increasing  Reynolds 
number. 
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Fig.  11  Shear-layer  mixed-fluid  fraction,  Sm/^T,  (local)  Reynolds  numbers.  Res, 
for  incompressible  flows  (Me  <  0.25).  Diamonds:  chemically-reacting  flows, 
squares:  cold-chemistry  flows.  Model  predictions:  solid  curve  -  BBM, 
dashed  curve  -  D87. 


3.  Mixing  and  the  geometry  of  isoscalars  in  turbulent  flows 


During  this  reporting  period,  an  analysis  was  undertaken  to  quantify  the  ge¬ 
ometry  of  isoscalar  surfaces  derived  from  scalar-field  data  in  turbulent  fiows.  In 
particular,  the  data  used  were  the  measurements  by  Catrakis  &  Dimotakis  (1996a) 
of  mixing  in  liquid-phase  (high-Schmidt-number)  turbulent  jets  and,  more  recently, 
three-dimensional  DNS  data  of  mixing  in  Rayleigh-Taylor-instability-driven  flows 
(Cook  1998). 

For  the  jet  data,  the  fluid  Schmidt  nximber  was  Sc  ~  2.0  x  10®,  with  flow 
Reynolds  numbers  of  Re  ~  4.5  x  10®,  9.0  x  10®,  and  18  x  10®.  These  fluid/flow 
parameters  resulted  in  a  jet-fluid  concentration  field  with  a  complex,  multiscale 
geometry.  Two-dimensional,  high-resolution  and  high-signal-to-noise-ratio  images 
of  the  scalar  far  field  in  this  flow  were  recorded  using  the  laser-induced-fluorescence 
technique,  in  a  plane  normal  to  the  jet  axis,  and  are  shown  in  Fig.  12.  They  span 
the  whole  jet-fluid-concentration  field,  at  the  z/dj  =  275  downstream  measurement 
location. 


Fig.  12  Laser-induced  fluorescence  scalar-field  data  in  a  liquid-phase  turbulent  jet 
at  Re  ~  4.5  x  10®  (left)  and  18  x  10®  (right)  in  a  plane  normal  to  the  jet 
axis  in  the  far  field  (z/d;  =  275). 
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It  is  interesting  to  note  that  classical  scalar  spectra  convey  little  information 
regarding  geometry  of  isoscalar  surfaces.  Spatial  spectra  computed  for  the  jet  scalar- 
field  data  are  shown  in  Fig.  13.  The  2-D  spectrum  (left)  is  for  a  single-image  realiza¬ 
tion  (Re  =  9  X 10^)  and  can  be  seen  to  be  very- nearly  axisymmetric.  Radial  spectra, 
obtained  by  azimuthal  integration,  are  also  shown  (right)  for  the  three  Reynolds 
numbers  investigated.  As  can  be  seen,  other  than  axisymmetry,  very  little  informa¬ 
tion  is  conveyed  by  the  2-D  spectrum.  A  similar  lack  of  information  can  be  noted 
for  the  radially-averaged  spectra,  other  than  a  decreasing  wavenumber  content  with 
increasing  Re,  in  contrast  to  classically-expected  behavior.  The  departures  in  be¬ 
havior  occurs  at  wavenumbers  corresponding  to  scales  roughly  1/3  the  image  extent 
and  are  not  the  consequence  of  measurement  resolution,  which  is  adequate. 


Fig.  13  Scalar  spectrum  for  single  realization  in  a  turbulent  jet,  at  Re  =  9  x  10® 
(left)  and  ensemble-averaged  radial  scalar  spectra  (right),  for  Re  =  4.5  x 
10®,  9  X  10®,  and  18  x  10®;  fines  of  increasing  solidity  denote  increasing  Re. 


For  the  jet  isoscalar-surface  analysis,  the  concentration  data  were  thresholded 
at  three  scalar  levels,  c(x,y)  =  cj,  C2,  and  C3,  where  the  C2  level  was  chosen  to 
correspond  to  the  peak  in  the  scalar  pdf  at  the  two  lower  Reynolds  numbers  and  the 
Cl  and  C3  levels  span  the  C2  level  and  represent  outer  and  inner  contours,  respectively 
(cf.  Catrakis  &  Dimotakis  1996a,  Fig.  8)  The  majority  of  the  level-set  analysis 
was  performed  at  the  C2  threshold.  One  of  the  most  interesting  findings  is  that 
the  coverage  length  of  the  isoscalar  contours,  defined  as,  L(X)/6b  =  AA^2(-^)/^b5 
decreases  with  increasing  Re.  This  is  plotted  in  Fig.  14  (left).  The  size,  ^b)  of  the 
bounding  box  partitioned  to  obtain  the  coverage  coxmts  was  independent  of  Re, 
within  measurement  statistics,  at  the  C2  scalar  threshold. 
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Another  finding  is  that  the  coverage  dimension,  D2(A),  is  a  smoothly-increasing 
function  of  scale,  boxmded  by  its  limiting  value  of  unity  (topological  dimension),  at 
the  smallest  scales,  and  2  (embedding  dimension),  at  the  largest  scales,  for  these 
data.  This  behavior  persists  for  all  Reynolds  numbers  investigated  and  is  important 
implying  that  previously-reported  power-law  coverage  proposals  seriously  misrep¬ 
resent  the  total  length  of  such  contours.  Accurate  knowledge  of  such  quantities  is 
crucial  for  the  prediction  and  modeUng  of  mixing,  chemical  reactions,  and  combus¬ 
tion  in  turbulent  flows. 

Both  the  coverage  length  and  dimension  indicate  less-convoluted  level  sets 
with  increasing  Re,  in  accord  with  the  radially- averaged  spatial-spectrum  sequence 
(Fig.  13,  right).  The  limiting  value  of  I'(A),  as  A  — »•  0,  and  the  coverage  dimension, 
D2(A),  at  meditnn-to-laxge  scales,  both  decrease  with  increasing  Re.  These  findings 
are  consistent  with  enhanced  mixing,  relative  to  stirring,  as  Re  increases,  leading  to 
improved  local  homogenization  of  the  scalar  field  and  geometrically-simpler  scalar 
level  sets.  These,  in  tmn,  result  in  lower  surface-to-volume  ratios,  with  increasing 
Re.  This  is  manifest  in  the  comparison  plot  of  sample  C2-level  sets  computed  from 
Re  =  4.5  X  10^  (left)  and  Re  =  18  x  10®  (right)  realizations  (Fig.  15). 


log,o  (Vi5b)  log,o  (V<5b) 


Fig.  14  Coverage  length  and  coverage  dimension  as  a  function  of  (normalized)  scale 
for  scalar  level  sets  in  a  turbulent  jet.  Re  ~  4.5  x  10®:  dotted/crosses; 
Re  ~  9.0  X  10®:  dashed/triangles;  Re  ~  18  x  10®:  solid/squares. 


A  new  framework  was  also  developed  for  the  study  of  the  geometry  of  level 
sets,  in  general,  that  extends  standard  (self-similar)  fractal  analysis.  Relations  be- 
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Fig.  15  Scalar  level  sets  at  Re  ~  4.5  x  10^  (left)  and  18  x  10^  (right),  indicating 
simpler  topology  at  higher  Reynolds  number. 

tween  coverage  statistics  and  scale  distributions,  in  general,  were  derived  (Catrakis 
&  Dimotakis  1996b).  The  scale  distribution  is  formulated  in  terms  of  a  multidimen¬ 
sional  (2-D,  for  these  data)  spacing  scale,  defined  as  the  size  of  the  largest  empty 
box  (LEB),  randomly  located  in  the  boimding  box  of  the  level  set,  that  contains 
no  part  of  the  level  set.  At  the  inner  scales,  the  coverage  dimension  and  LEB- 
scale  distribution  were  found  to  be  consistent  with  lognormal  statistics.  The  LEB 
scale  distribution  was  also  shown  to  provide  a  dimensionless  measure  of  the  surface- 
volume  (perimeter-area,  in  2-D)  ratio  of  the  level  sets,  in  the  small-scale  limit.  A 
more  detailed  description  of  scale  distributions  and  their  relation  to  (power-law  or 
scale-dependent)  fractal  geometry  can  be  found  in  Catrakis  &  Dimotakis  (1996c). 

For  the  jet  data,  level  sets  consist  of  individual  (disjoint)  “islands”  and  “lakes”, 
depending  on  whether  the  interior  is  at  a  lower,  or  higher,  scalar  level,  respectively. 
It  is  useful  to  analyze  island/lake  statistics,  such  as  size  and  shape  complexity.  In  the 
context  of  combustion,  for  example,  an  island  would  be  associated  with  an  unbumt 
fuel  pocket  in  a  non-premixed  turbulent-jet  fiame.  Such  an  analysis  indicates  that 
the  size  distribution  of  such  features  is  well  approximated  by  a  log-normal  pdf,  at 
small-to-intermediate  scales  (Catrakis  Sz  Dimotakis  1996b),  as  shown  in  Fig.  16. 
Size  here  is  defined  as  \/A,  with  A  the  individual  island/lake  area. 

Returning  to  chemical  reactions  and  combustion  in  non-premixed  hydrocarbon 
turbulent  flames,  in  which  combustion  is  largely  confined  to  the  instantaneous  sto- 
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ichiometric  (isoscalar)  surface  (Burke  &  Schumann  1928),  area-volume  measures  of 
the  isoscalar  surface  may  be  used  to  relate  the  local  burning  rate  to  the  time  re¬ 
quired  for  the  local  consumption  of  unbumt  fuel  pockets.  Such  a  measure,  dubbed 
shape  complexity,  can  be  defined  as,  1  <  ^2  =  ■P/[2(7r.4)^/^],  in  2-D,  where  P  is 
the  perimeter  and  A  the  area  of  an  island  or  lake,  with  (^2)jciin  ~  ^  attained  for 
a  circle,  and  corresponding  extensions  for  0,^,  for  higher-dimensional  embedding 
spaces.  The  liquid-phase  jet  data  described  above  indicate  that  a  power  law  over 
3  decades  in  size  (6  decades  in  area)  provides  a  good  approximation  for  the  pdf 
of  shape  complexity.  This  behavior  is  equivalent  to  log-Poisson  statistics  for  ^2 
(Catrakis  &  Dimotalcis  1998). 

A  coverage  analysis  was  also  performed  on  isodensity  data  from  a  Navier-Stokes 
DNS  study  of  the  evolution  of  a  Rayleigh- Taylor-instability  flow,  of  a  Sc  =  u/V  =  oo 
fluid.  The  flow  was  initialized  with  a  p  =  3  fluid  on  top  and  a  p  =  1  fluid  on  the 
bottom,  in  a  256  x  256  x  512  rectangular  box.  The  three-dimensional  DNS  of 
the  evolving  flow  was  terminated  when  the  spatial-resolution  requirements  could 
no  longer  be  met  by  the  fixed  grid,  at  Rcfinal  —  1-1  x  10^,  based  on  the  vertical 
extent  and  growth  rate  of  the  Rayleigh-Taylor  mixing  region  (Cook  1998).  The 
simulation  utilized  periodic  boimdary  conditions  in  the  boundary  planes  transverse 
to  the  acceleration  vector,  and  no-slip  at  the  top  and  bottom  faces  at  the  end  of  the 
long  dimension  of  the  box;  aligned  with  the  acceleration  vector.  A  small- amplitude 
perturbation  of  the  interface  between  the  two  fluids  initialized  the  flow. 
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Fig.  17  Temporal  evolution  of  scalar  power  spectrum  (left)  and  coverage  dimension 
(right)  for  a  2D  slice  (at  midheight)  of  the  p  =  2  isoscalar  surface  in  a 
numerical  simulation  of  the  Rayleigh- Taylor-inst  ability  flow. 

As  with  the  jet-flow  data,  spectral  information  is  not  enough  to  characterize 
the  geometry  of  isoscalar  surfaces.  Figure  17  (left)  shows  the  time  evolution  of  the 
spatial  spectrum  of  the  density  field  in  the  midheight  plane,  i.e.,  of  p{x,y,z  =  0), 
for  the  indicated  progression  in  time.  The  spectrum  initially  exhibits  a  tempo¬ 
ral  progression  to  lower  wavenumbers,  a  result  of  diffusive  smoothing  of  the  initial 
density-field  perturbation  (recall  that  Sc  =  1  here).  While  the  low  wavenumber 
spectral  content  continues  to  increase,  a  sustained  progression  to  higher  wavenum¬ 
bers  with  increasing  time  (for  t  >  1.4),  can  be  seen,  plausibly  as  the  Rayleigh- Taylor 
and  other,  secondary,  flow  instabilities  take  over,  with  the  spectrum  reflecting  the 
growth  of  small-scale  features  at  the  midheight  plane. 

Figure  17  (right)  plots  the  coverage  dimension,  D2{X),  for  the  p(x,  y,  z  =  0)  =  2 
isodensity  contomrs,  in  the  same  plane,  computed  by  successive  binary  subdivisions 
of  the  midheight  slice,  as  in  the  jet  scalar-data  analysis,  above.  The  resulting 
scale- dependent  coverage  dimension  D2{X)  can  be  seen  to  span  the  range  of  values 
from  unity  (the  topological  dimension),  to  2  (the  embedding  dimension),  smoothly 
transitioning  between  the  two  limiting  values,  at  the  smallest  and  the  largest  scales, 
respectively.  Interestingly,  the  temporal  progression  indicated  by  D2{y)  is  from 
small  to  large  scales,  for  t  <  3.5,  i.e.,  opposite  the  high-wavenumber  trend  in  the 
spectral  analysis.  It  is  not  until  near  the  end  of  the  simulation  (for  t  >  3.5),  where 
a  reversal  of  this  trend  is  exhibited,  at  small  scales  only. 
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The  spectnim  measures  the  wavenumber  content  of  the  selected  density  field, 
p{x,y,z  =  0;<),  i.e.,  of  the  density  surface  over  the  midheight  (x,  y)-plane,  while 
the  level-set  coverage  analysis  interrogates  the  geometry  of  the  p(x,y,z  =  0;t)  =  2 
contours  in  the  midheight  (x,y)-plane,  i.e.,  the  slice  of  the  former,  at  the  p  =  2 
elevation.  The  early-time  (t  —  0.938  and  1.875)  contomrs  become  smoother,  while 
the  density  surface  is  developing  sharp  peaks.  The  latter  are  responsible  for  the 
spectral  progression  to  higher  wavenumbers.  This  trend  continues  for  a  while,  with 
a  subsequent  transition  that  can  be  seen  in  the  behavior  of  the  density  surfaces. 
These  develop  a  more  complicated  topology,  characterized  by  folds  between  local 
maxima  and  minima,  at  late  times  (t  —  4.607  and  4.810).  This  also  registers  in  the 
level-set  contours  that  can  be  seen  to  develop  small-scale  features,  on  top  of  the 
larger-scale  features  that  continue  to  increase  in  size. 

Viewing  the  flow  evolution  through  the  three-dimensional  density-field  data 
(not  discussed  here)  indicates  that  the  likely  cause  of  this  transition  is  the  devel¬ 
opment  of  secondary  instabilities,  of  the  Kelvin-Helmholtz  type  in  the  high-shear, 
near-midheight  regions  generated  by  the  interpenetrating  Rayleigh- Taylor  fingers, 
and  the  formation  of  mushroom-like  structures  at  their  tips.  The  cross-over  in  the 
coverage  dimension,  Z)2(A),  at  small  scales  for  late  times  revealed  this  transition, 
even  though  there  is  scant,  if  any,  evidence  for  it  in  the  spectral  data.  It  is  an 
attestation  of  the  scale-local  capability  of  the  coverage  analysis  that  such  geomet¬ 
rical  properties  were  clearly  registered  in  those  statistics.  They  were  subsequently 
confirmed  by  computer-visualization  of  the  corresponding  field  information. 

Additional  details  of  this  part  of  the  work  were  presented  at  the  2“*^  Monte 
Verita  Colloquium  on  Fundamental  Problematic  Issues  in  Turbulence,  22-28  March 
1998  (Ascona,  Switzerland)  and  documented  in  Dimotakis  et  al.  (1998a). 


4.  Image  correlation  velocimetry 

The  brief  account  below  is  abstracted  from  the  thesis  by  Gornowicz  (1997), 
which  should  be  consulted  for  additional  documentation. 

The  Image  Correlation  Velocimetry  (ICV)  method  developed  as  part  of  work 
sponsored  by  this  grant,  is  an  extension  of  the  original  implementation  by  Tokumaru 
&  Dimotakis  (1995).  The  new  implementation  utilizes  a  hierarchical  B-spline  rep¬ 
resentation  of  the  velocity  (mapping)  field,  allowing  the  desired  spatial-continuity 
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order  to  be  imposed  on  the  solution,  as  will  be  discussed  below.  A  very  simi¬ 
lar  approach  is  used  by  Szeliski  &:  Shum  (1996).  These  authors  also  employ  a 
multiresolution-spline  representation  of  the  displacement  field  between  image  pairs. 

The  development  of  the  ICV  extensions  described  here,  is  very  similar  to  the 
Szelinski  &  Shum  method,  that  was  developed  in  a  different  context.  Variational 
methods  have  been  used  to  infer  displacement  fields  in  a  variety  of  contexts,  as 
done  by  Zhou  et  ah  (1995),  for  example,  who  employed  a  multiresolution  repre¬ 
sentation  of  the  three-dimensional  displacement  field  in  the  interior  of  a  cylindrical 
asphalt /aggregate  core,  assuming  a  volume-preserving  (divergence-free)  displace¬ 
ment  field. 

Another  general  category  of  this  methodology  identified  by  Barron  et  al.  (1994) 
are  the  so-called  differential  techniques,  pioneered  by  Horn  &  Schimck  (1981).  These 
methods  calculate  the  components  of  the  scalar-transport  equation  and  use  addi¬ 
tional  constraints  to  remove  ambiguities.  Recent  work  has  investigated  and  com¬ 
pared  the  required  additional  constraint(s)  to  the  scalar  transport  equation  pro¬ 
posed  by  various  researchers  after  Horn  &  Schxmck  (e.p.,  Willick  k.  Yang  1991). 
Strong  proponents  of  the  application  of  this  technique  for  fluid  velocimetry  have 
been  Dahm  et  al.  (1991,  1992).  More  recently,  a  variational  approach  was  offered 
by  Su  k  Dahm  (1995)  and  Dahm  et  al.  (1996).  Pearlstein  k  Carpenter  (1995), 
however,  noted  that  the  method  of  Dahm  and  collaborators  suffers  from  a  local 
ambiguity  problem  in  that  the  local  velocity  field  is  only  defined  in  the  direction  of 
the  imaged-scalar  gradient.  Pearlstein  k  Carpenter  proposed  to  mitigate  this  am¬ 
biguity  problem  by  simultaneously  tracking  multiple  scalar  fields.  A  shortcoming  of 
these  methods  is  that  they  rely  on  spatial  and  temporal  derivatives  of  the  imaged 
field  to  deduce  the  convecting  velocity  field,  rendering  them  rather  susceptible  to 
the  inevitable  image  noise. 


4.1  Continuous-field  ICV  method 

The  ICV  procedure  seeks  the  displacement  field,  ^(x),  such  that  the  region 
in  the  neighborhood  of  x,  in  the  image  Ii(x),  at  time  ti,  is  best  mapped  into  the 
region  x  -f  ^  in  the  next  image,  i2(x),  at  t2  —t\-\-  r,  i.e.,  such  that, 


Ji(x)  1-^  /2(x  +  0  . 


(7) 
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Specifically,  the  displacement  field,  ^(x),  is  sought  that  minimizes  the  square  of  the 
difference  of  the  two  images,  integrated  over  the  correlation  domain,  fi,  z.e.,  a  cost 
function  given  by, 

=  j [hi'x.  -  Ii{-x.)f  df2(x)  min  .  (8) 

n 

This  cost  fimction  alone  does  not  guarantee  either  a  unique  or  smooth  solution. 
Such  attributes  depend  on  the  functional  representation  of  ^  and  are  addressed  in 
Gomowicz  (1997). 

If  the  time  difference,  r  =  <2  -  j  between  the  two  images  is  small,  in  some 
appropriate  sense,  one  can  Taylor-expand  the  displaced-image  field  at  <2,  *-e., 

72(x  +  ^,ti+T)  =  7i(x,ii)  +  T^Ii(x,ti)  +  ^- ^/i(x,ti)  +  H.O.T.'s  . 

where  the  higher-order  terms  are  O  (r^),  O  (^^),  or  A  mapping  (displace¬ 

ment)  field  (Eq.  7),  i.e.,  one  that  produces, 

/2[x-|-^(x),<i -f-r]  ~  /i(x,ti)  , 
is  seen  to  be  equivalent  to  the  requirement  that, 

again,  to  leading  order  in  the  space/time  displacements,  i.e., 

d  d  1 

^  ^  ~  ^  ^  for,  u  =  -  ^  .  (9a,b) 

Equation  9  is  the  optical-flow  equation.  It  is  also  the  scalar-transport  equation, 
provided  diffusive  effects  are  negligible,  which  typically  translates  to  an  upper  limit 
on  the  time  interval,  r ,  between  the  image  pair.  Since  scalar  diffusivity  is  essentially 
fixed  by  the  choice  of  the  fluid,  the  time  interval  must  be  chosen  such  that  diffusion 
is  neghgible  (cf.  discussion  in  Tokumaru  &  Dimotakis  1995). 

To  compute  an  optimal  mapping  field,  the  ICV  method  rehes  on  a  parametric 
representation  of  the  displacement  field,  ^(x).  In  several  refined  DPIV  implemen¬ 
tations  (e.g.,  Huang  1994,  Sholl  &  Savas  1997),  as  well  as  in  the  previous  ICV 
implementation  (Tokumaru  &  Dimotakis  1995),  local  Taylor  expansions  of  the  dis¬ 
placement  field  were  employed,  to  various  orders,  i.e., 

$(x)  =  ^(xc)  +  (x-Xc)iac,i 

+  I  (X  -  Xc)i  (X  -  Xc)j /3c,ij 

"b  ^  Xc)i  (x  Xc)j  (x  Xc)^  'fcjijk 
+  etc.  , 


(10) 
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where  Xc  is  any  expansion  point,  with  the  parameters,  ac,i,  ^c,iji  7c,ijk,  deter¬ 
mined  by  the  DPIV  or  ICV  solution. 

In  the  ICV  implementation  of  Tokumaxu  &  Dimotakis  (1995),  the  minimized 
cost  function  included  terms  which  increased  its  value  with  the  (square  of  the) 
amplitude  of  any  discontinuities  of  the  displacement  field  and  its  derivatives  at 
the  bovuidaxies  of  the  array  of  Taylor-expansion  regions  around  the  selected  con¬ 
trol  points,  Xc-  As  a  consequence,  much  of  the  built-in  flexibility  in  describing 
spatial  variations  of  the  displacement  (velocity)  field  was  lost,  with  degrees  of  free¬ 
dom  gained  from  the  Taylor-expansion  coefficients  in  Eq.  10,  in  effect,  expended  to 
minimize  discontinuities  of  the  velocity  field  and  of  its  derivatives  at  the  Taylor- 
expansion  region  boundaries. 

To  mitigate  this  difficTilty,  the  present  ICV  implementation  relies  on  a  displace¬ 
ment  field  that  possesses  the  required,  C”,  continuity  properties  by  construction, 
with  the  order  of  continuity,  n,  chosen  appropriately,  as  described  below.  The  re¬ 
maining  (true)  degrees  of  freedom  axe  utilized  to  minimize  the  cost  function, 
with  no  added  (smoothing)  terms  in  the  integrand. 

Velocity-  and  vorticity-field  solutions  of  the  Navier-Stokes  equations  are  con¬ 
tinuous,  with  continuous  derivatives  to  all  orders,  i.e.,  axe  C°°.  In  the  present 
implementation,  which  was  limited  to  two-dimensional  fields,  a  displacement 
(velocity)  field  was  employed,  i.e.,  possessing  continuous  second  derivatives,  cor¬ 
responding  to  inferred  vorticity  fields  that  possessed  continuous  first  derivatives. 
This  was  achieved  by  representing  the  displacement  field  in  terms  of  B-sphnes  with 
appropriate  basis  fimctions,  whose  control  parameters,  G  3?^,  then  provided 

the  parametric  description  of  the  displacement  field,  i.e., 

«x)  =  {  [x;  qS;’''*]  .  (11) 


With  the  solution  space  of  the  minimization  problem  (Eq.  8)  restricted  in  this 
fashion,  the  cost  functional,  becomes  a  function  of  the  control  parameters. 


I.e., 


J 


Ar,R) 

ilj 


(12) 


(r  K) 

possessing  a  minimxun  where  derivatives  of  ff,  with  respect  to  each  of  the  q)j’  , 
vanish.  This  allows  a  global  minimization  over  the  (selected)  image-correlation  do¬ 
main  to  be  sought,  using  an  iterative,  multidimensional,  conjugate-gradient  method 
over  all  parameter  values,  with  a  suitable  initial  guess. 
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4.2  ICV  algorithm  implementation 

The  ICV  implementation  is  comprised  of  a  sequence  of  iterative,  algorithmic 
steps: 

•  image-data  preparation, 

•  image-correlation  domain  definition, 

•  cross-correlation  displacement-field  initialization,  and 

•  a  final,  conjugate-gradient  displacement-field  optimization  to  minimize  the 
cost  function  (Eq.  8). 


The  procedure  starts  by  further  processing  individual  data  images,  after  back¬ 
ground  removal,  illumination  normalization,  etc.,  for  shot-to-shot  intensity  varia¬ 
tions  of  the  illuminating  laser  sheet.  A  geometry  file  is  generated,  which  specifies 
the  correlation  domain,  fi,  within  the  image  domain.  An  initial  hierarchy  of  the  B- 
spline  resolution  knot  grid  is  specified  and  any  excluded  regions  from  the  correlation 
domain  {e.g.,  laser  shadows,  imaging  occlusions,  etc.)  are  also  identified.  The  outer 
botmdary  of  the  correlation  domain,  fi,  is  specified  as  a  polyline  (n-sided  polygon). 
Inner  boundaries  can  also  be  accommodated,  allowing  correlation  domains  to  be 
defined  that  axe  not  simply-connected,  as  necessary.  These  boundaries  are  defined 
on  the  first  of  the  two  images  (for  each  pair). 

The  next  step  initializes  the  solution  at  the  coarsest  resolution  level;  usually, 
one  spline  patch.  The  initialization  is  performed  by  cross-correlating  spatially-local 
windows,  using  Fourier  techniques,  as  in  DPIV  analyses  (e.^.,  Adrian  1991,  Willert 
&  Gharib  1991).  The  results  of  these  correlations  initialize  the  mapping  vector 
field,  ^(x).  Windows  are  then  centered  at  the  peak  of  each  B-spline  basis  function 
and  the  cross-correlation  results  are  used  to  determine  the  corresponding  B-spline 
control  parameters,  at  the  resolution  level  r. 

The  initialization,  of  the  B-spline  representation  for  allows  Eq.  7  to  be 
applied,  producing  an  initial  mapped  version  of  the  second  image,  i.e.,  I2{x  + 
that  is  “closer”  to  the  first  image,  Ii(x).  Further  cross-correlations  are  run  between 
/i(x)  and  /2(x  -1-  ^)  to  produce  subsequent  estimates, 

A  similar  process  for  determining  the  displacement  field  in  DPIV,  but  without 
FFT’s,  is  outlined  in  Huang  (1994)  and  termed,  “Particle  Image  Distortion”.  A  fast 


25 


version,  termed,  “Lagrangian  Particle  Tracking”,  was  introduced  by  Sholl  &  Savas 
(1997).  In  these  implementations,  the  deduced  displacement  field  was  specified  in 
terms  of  local  Taylor  expansions,  to  first  and  second  order,  respectively. 

The  ICV  initialization  step  starts  with  large  cross-correlation  windows  (up  to 
256  X  256  pixels,  or  1/4  x  1/4  the  image  domain,  for  example)  to  avoid  spurious 
correlations  and  to  pick  up  any  large  displacements.  This  initialization  step  is  par¬ 
ticularly  important,  if  there  are  displacements  greater  than  half  the  characteristic 
length  of  a  continuous  scalar  used  to  mark  the  flow  (equivalent  to  a  Nyquist  crite¬ 
rion).  Small-scale  features  of  the  velocity  field  are  determined  in  subsequent  stages. 
This  aspect  is  particularly  important,  in  as  much  as  the  subsequent  minimization 
stages  may  not  correct  for  errors  introduced  at  this  stage  and  a  local  minimum  of 
J (^)  might  be  foimd  instead. 

Once  large-scale  displacements  have  been  found  with  such  windows,  the  size 
of  the  window  is  successively  reduced  by  a  factor  of  2,  cross-correlations  are  per¬ 
formed,  and  the  corresponding  B-spline  parameters  are  computed  to  yield  the  next 
window-size  estimates  of  the  displacement  (mapping)  field.  These  successive  halv¬ 
ings  continue  until  a  user-determined  minimum  window  size  is  reached. 

The  displacement  field,  produced  by  the  cross-correlation  sequence  is  used  to 
initialize  an  iterative  minimization  procedure.  The  projection  of  the  displacement 
field  on  the  set  of  B-spline  basis  functions  (c/.  Eq.  11  and  discussion  in  Gornowicz 
1997)  converts  the  integral  J  to  be  minimized  from  a  functional  of  ^  to  a  function 
of  the  finite  number  of  B-spline  control  parameters  (c/.  Eq.  12),  as  noted  above. 
The  required  minimization  of  J  is  now  performed  in  a  finite- dimensional  space. 

An  important  part  of  the  ICV  implementation  is  the  accommodation  of  ap¬ 
propriate  boundary  conditions.  The  inferred  flow  fields  are  typically  elliptic,  so  the 
choice  of  appropriate  boundary  conditions  is  important.  Especially  important  is 
the  inference  of  flows  in  the  vicinity  of  no-slip  boundaries,  an  area  that  is  typically 
not  very  well  resolved  by  DPIV-correlation  methods. 
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5.  Three-dimensional  imaging  and  visualization 

As  part  of  an  effort  to  investigate  the  three-dimensional  structure  of  scalar 
fields,  a  set  of  three-dimensional,  laser-induced  fluorescence  (LIF)  measurements  of 
scalar-species  concentration  in  the  far-field  of  a  liquid-phase  axisymmetric  turbulent 
jet  was  undertaken.  The  spatio-temporal  data,  consisting  of  time  sequences  of  LIF 
images,  capture  the  entire  local  extent  of  the  jet  and  are  spatially  resolved  to  the 
diffusion  length  scale  at  Re  ~  2.3  x  10^.  An  in-house  designed  and  constructed 
digital  camera  was  utilized  to  record  high-resolution  (1024  x  1024),  high-dynamic- 
range  (12-bits/pixel)  images,  acquired  at  a  rate  of  lOframes/s. 

The  space-time  structure  of  the  turbulent  mixing  flow  was  visualized  by  recon¬ 
structing  three-dimensional  isoconcentration  surfaces  from  the  sequence  of  LIF  im¬ 
ages.  In  collaboration  with  Dr.  David  Laidlaw,  of  the  Caltech  Biology  and  Computer 
Science  Department,  isosurfaces  for  a  390  x  390  x  4,  3-D  portion  of  the  turbulent 
jet  data  were  calculated  and  visualized  using  the  marching-cubes  algorithm  com¬ 
monly  used  in  medical  imaging  (e.g.,  computed  tomography,  magnetic  resonance, 
and  single-photon-emission  computed  tomography). 

For  high  Reynolds  number  turbulence,  it  is  all  but  impossible  to  satisfy  strict 
pixel-by-pixel  application  of  the  Nyquist  sampling  criterion.  As  a  result,  the  con¬ 
ventionally-constructed  isoconcentration  surface  in  Fig.  18  is  rough,  exhibiting  steps 
on  the  front  surface  as  well  as  cone-like  structures  along  the  left  side  of  the  image 
that  are  artifacts  of  temporal  aliasing.  This  aliasing  occurs  in  the  standard  3-D 
isosurface  visualization  despite  the  fact  that  the  LIF-image  data  are  temporally 
resolved  with  respect  to  the  diffusion-scale  advection  time  —  Aj)/I7. 

Image  Correlation  Velocimetry  (ICV),  as  described  above,  extracts  velocity 
information  from  a  succession  of  images  of  convected  scalars  and  can  be  applied 
to  advantage  to  LIF  scalar-field  data  to  facilitate  three-dimensional  visualization. 
ICV-derived  velocity  information  is  used  to  generate  intermediate  interpolated  (su- 
persampled)  2-D  scalar  images  that  can  then  be  used  by  standard  visualization 
algorithms  to  reconstruct  the  three-dimensional  structure  of  scalar  fields. 

To  satisfy  the  Nyquist  criterion  for  temporal  sampling,  ICV  was  used  to  con¬ 
struct  an  additional  19  images  in-between  each  successive  frame  of  LIF  data.  As 
discussed  above,  the  displacement  field,  ^(x),  is  represented  using  bicubic  B-splines 
at  a  hierarchy  of  spatial  resolutions  (i.e.,  multiresolution)  and  the  minimization  of 
J7{^(x)}  (Eq.  8)  is  implemented  with  a  multidimensional  conjugate-gradient-based 


Fig.  18  3-D  space-time  rendering  of  scalar  isoconcentration  surface  derived  from 
LIF  measurements  in  the  far-field  of  a  tiurbulent  jet  at  Re  ~  2.3  x  10* 

approach.  With  this  information,  temporally-interpolated  images  are  computed 
based  on  the  ICV-calculated  displacement  field,  ^(x),  as, 

(x)  =  ^  { -fi  (x  -  +  h  [x  -  (1  -  a) ^]}  , 

where  i  is  an  integer,  1  <  i  <  19,  and  a  =  i/20. 

Figure  19  depicts  an  isoconcentration  surface  computed  with  the  marching- 
cubes  algorithm  from  the  same  LIF  data  as  in  Fig.  18,  using  the  ICV-interpolated 
intermediate  images.  The  dealiased  (20  x  supersampled)  three-dimensional  space- 
time  isosurface  is  noticeably  smoother  than  the  original  visualization  (Fig.  18).  An 
inclined  cylindrical  isosurface  can  be  seen  along  the  left  edge  of  the  image,  rather 
than  two  pairs  of  disconnected  “cones”.  Evidence  of  developing  Kelvin-Helmholtz 
instabilities  can  now  be  seen  on  the  right  half  of  the  isoconcentration  surface  in 


Fig.  19  Dealiased  3-D  rendering  of  scalar  isoconcentration  surface  (same  data  as 
in  Fig.  18)  using  Image  Correlation  Velocimetry  to  generate  intermediate 
images  in  time. 

Fig.  19.  The  inverse-slope  and  orientation  of  the  structure  correspond  to  the  local 
convection  velocity  of  the  isoconcentration  surface.  Our  recently-developed  tech¬ 
nique  of  incorporating  ICV  into  the  isosurface  computation  has  effectively  increased 
the  temporal  resolution  of  LIF  images  of  the  scalar  field. 


6.  Chemical-kinetic  simulations 

As  described  above,  the  experimental  investigations  of  molecular  mixing  in  sub¬ 
sonic  and  supersonic  turbulent  shear  layers  have  relied  on  (H2  H-N0)/F2  hypergolic 
mixtures  in  passive  diluents  (N2,  Ar,  He,  etc.).  This  chemical  system  is  hypergolic, 
with  H2  and  NO  premixed.  On  contact,  NO  and  F2  react  spontaneously,  producing 
F  radicals,  that  break  the  H2  molecules. 
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The  (H2  +  N0)/F2  system  can  be  kinetically-fast  even  at  modest  reactant 
concentrations  and  will  react  to  produce  a  temperature  rise  that  can  be  directly 
related  to  the  number  of  moles  of  chemical  product  in  the  timbulent  shear  layer 
(Dimotakis  1991a).  In  this,  high-Damkohler-number  limit,  the  rate  of  chemical 
reaction  and  heat  release  is  limited  by  the  rate  of  molecular  mixing. 

Chemically-reacting  mixing-layer  experiments,  such  as  those  described  above, 
indicate  that  the  normalized  temperature  rise, 

^  T-Too 

~  Ti-T^  ’ 

where  Tf  is  the  adiabatic  flame  temperature  and  Too  the  temperature  of  the  imre- 
acted  mixture,  reaches  peak  (mean)  values  in  mixing  layers  in  the  range, 

0.45  <  (^)„,ax  <  0.6  ,  (14) 

despite  the  fact  that  a  fraction  of  the  fluid  is  comprised  of  (as  yet)  unmixed  re¬ 
actant  freestream  fluids  (Dimotakis  1991a)..  Additionally,  the  relatively-low  con¬ 
centrations  of  the  (H2  -l-  N0)/F2  system  reactants  used  in  these  experiments  (and 
the  associated  simulations  described  below)  resulted  in  mixing-limited  (high-Da) 
product  formation. 

An  attempt  to  model  and  assess  the  effect  of  the  detailed  chemical  kinetics  in 
such  experiments  was  made  by  Hall  &  Dimotakis  (1987)  by  considering  a  simple 
model  in  which  the  chemical  reactions  were  assumed  to  proceed  in  a  perfectly- 
stirred  reactor.  This  approach  omitted  the  effects  of  fluid  mechanics  (strain  rate) 
and  molecular  transport. 

The  present  investigation  was  undertaken  to  augment  the  previous  studies  by 
conducting  a  detailed  numerical  simulation  of  strained,  laminar,  opposed-jet  react¬ 
ing  flows  for  the  (H2  -f-  N0)/F2  system,  to  assess  some  of  the  mechanisms  affecting 
the  details  of  the  structure  and  dynanoics  of  the  resulting  reaction  fronts.  This  flow 
configuration  provides  a  simple  description  of  the  behavior  of  local  diffusion/reaction 
surfaces,  in  the  turbulent  shear  zone,  that  retains  the  leading-order  features  of  the 
turbulence- chemistry  interactions. 

Numerical  solutions  for  a  counterflow  configuration  were  obtained  using  well- 
established  formulations  and  computer  codes  (Kee  et  al.  1988,  Kee  1991,  Kee  et 
al.  1985,  and  Kee  et  al.  1983)  .  The  counterflow  code  (Kee  1991)  was  further 
modified  to  accommodate  non-premixed  and  premixed  reacting  modes. 
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Aerodynamic  strain  rate, a,  s‘^ 


Fig.  20  Effect  of  aerodynamic  strain  on  maximum  temperature,  Tmax,  for  the  non- 
premixed  configuration. 

Non-premixed  cases  were  simulated  as  an  impinging  F2 /inert-jet  onto  an  oppos¬ 
ing  H2  +  NO/inert-jet  stream.  Premixed  combustion  modes  were  also  studied,  with 
identical,  opposing  (H2  +  NO)/F2 /inert  jets  impinging  on  each  other.  The  premixed 
mode  provides  an  approximate  description  of  the  combustion  of  premixed  reactants 
fed  to  the  reaction  zone  as  the  result  of  leakage  from  a  neighboring,  high-strain- 
rate  diffusion  layer;  a  likely  scenario  in  low-Damkohler  number,  turbulent-mixing 
flows.  The  kinetic  scheme  includes  nine  elementary  reactions  and  kinetic  parameters 
(Egolfopoulos  et  al.  1996). 


6.1  Non-Premixed  Mode 

Numerical  simulations  were  conducted  for  a  wide  range  of  strain  rates  and 
for  nine  different  reactant  cases,  defined  in  terms  of  reactant  concentrations  in 
their  respective  streams.  The  molar  concentrations  were  typical  of  those  used  in 
supersonic  mixing-layer  experimental  studies. 
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As  the  strain-rate,  cr,  increases,  the  maximum  temperature,  Tmax?  falls,  as  in 
high-activation  systems;  cr  was  defined  as  the  maximum  velocity  gradient  in  the 
convective  zone  of  the  F2  jet  stream.  However,  turning-point  extinction  does  not 
occur  in  such  mixtures  and,  while  Tmax  continues  to  decrease  with  increasing  cr ,  fi¬ 
nite  chemical  activity  persists  at  strain  rates  as  high  as  tens  of  thousands  reciprocal 
seconds.  This  is  illustrated  in  Fig.  20  for  a  Fa-rich  and  a  Fa-lean  reacting  config¬ 
uration,  in  terms  of  the  variation  of  the  normalized  temperature  rise,  9  (Eq.  13), 
with  cr.  In  these  simulations.  Too  was  300  K,  with  T{  computed  from  equilibrium 
calculations  using  the  STAN  JAN  program  (Reynolds  1986). 

The  mechanisms  responsible  for  the  temperature  reduction  with  increasing 
strain  rate  (c/.  Fig.  20)  were  also  studied.  Reaction-path  analyses  confirmed  that 
the  reaction, 

NO  Fa  =  NFO  +  F  (Rl) 

is  responsible  for  the  chain  initiation  by  releasing  an  F  atom.  This  is  used  by 
reaction, 

F  Ha  =  HF  -h  H  (R4) 

for  Ha  consumption  and  the  production  of  the  final  HF  product  and  active  H  radi¬ 
cals.  The  H  radicals  then  react  with  Fa,  through 

H  +  Fa  =  HF  +  F  ,  (R3) 

to  produce  HF  and  F,  which  is  further  used  in  R4. 

The  react ant-lealcage  mechanism  is  traceable  to  the  finite  kinetic  rate.  Char¬ 
acteristic  time  scales  for  the  destruction  of  NO,  Fa,  and  Ha  at  the  location  of  the 
maximum  rate  of  Rl,  R3,  and  R4,  respectively,  were  determined  and  scaled  by  the 
strain  rate,  <7,  to  form  Damkohler  numbers  for  all  three  species.  The  results  are 
shown  in  Fig.  21,  for  both  the  Fa-rich  and  Fa-lean  cases.  It  can  be  seen  that  Da’s 
decrease  with  a  and  attain  values  of  order  unity  at  high  strain  rates.  For  both 
Fa-rich  and  -lean  cases,  the  Da  of  Fa  have  lower  values  compared  to  NO  and  Ha. 
This  indicates  that  the  main  leakage  mechanism  is  that  of  Fa,  through  R3. 

These  observations  were  further  supported  by  detailed  sensitivity  analyses,  in 
which  the  effect  of  individual  elementary  reactions  on  Tmax  was  determined  at  dif¬ 
ferent  strain  rates.  Representative  results  are  shown  in  Fig.  22,  for  the  Fa-rich  and 
Fa-lean  cases.  For  both  cases,  Rl,  R3,  and  R4  favor  Tmax,  with  R3  having  the  high¬ 
est  positive  sensitivity.  Furthermore,  both  cases  confirm  that  Rl  is  very  important, 
initiating  the  reaction  between  the  two  main  reactants,  NO  and  Fa.  As  expected, 
negative  sensitivity  on  Tmax  is  exhibited  by  the  three-body  recombination  reactions. 


Damkohler  number  Damkohler  number 
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6.2  Premixed  Mode 

For  all  nine  non-premixed  configurations  and  for  strain  rates  of  the  order  of 
10^/s,  or  so,  which  axe  of  the  order  encountered  in  supersonic  mixing-layer  exper¬ 
iments,  Tmax  is  substantially  reduced  (Egolfopoulos  et  al.  1996).  This  indicates 
reduced  chemical  activity  and  that  unreacted  species  will  mix  and  react,  within 
highly-strained  flow  structures,  in  a  premixed  mode.  To  study  this  scenario,  var¬ 
ious  configurations  of  premixed  reaction  in  strained  flow  fields  were  investigated. 
For  those  studies,  the  composition  of  the  premixed  species  was  determined  using 
a  value  of  the  (molar)  entrainment  ratio,  between  the  high-speed  and  low-speed 
stream  (Dimotakis  1991b),  of  =  1.06  (Hall  et  al.l991  and  Bond  &  Dimotakis 
1993). 


1.00 
0.90 
H  0.80 

I 

•S  0.70 

H 

^  0.60 
0 

^  0.50 

I  0.40 
Y  0.30 
S  0.20 
0.10 
0.00 

10  100  1000  10000  100000 

Aerodynamic  (global)  strain  rate,  Ggiobai,  s'^ 

Fig.  23  Effect  of  aerodynamic  strain  rate  on  maximum  temperatme,  Tmax?  for  the 
premixed  configuration. 
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Two  identical  premixed  jets  impinging  onto  each  other  were  studied.  As  a 
consequence  of  this  implementation,  for  low  strain  rates,  losses  resulting  from  high 
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temperature  gradients  at  the  nozzle  exit  result  in  Tmax  well  below  Tf.  As  the 
strain  rate  increases,  temperature  gradients  at  the  nozzle  exit  decrease,  decreasing 
boundary  losses,  and  Tmax  approaches  Tf,  corresponding  to  homogeneously-mixed 
reactants  (cf.  Fig.  23).  In  the  premixed  mode  simulations,  the  convective  transport 
was  measured  by  a  global  strain  rate,  a^,  defined  as  the  nozzle-exit  velocity,  divided 
by  half  the  distance  between  the  two  nozzles.  For  all  other  cases  studied,  the 
response  of  Tmax  to  the  imposed  strain  rate  was  foxmd  to  be  similar. 

The  asymptotic  behavior  for  high  strain  rates,  i.e.,  Tmax  ~^T{,is  at  variance 
with  typical,  high-activation,  premixed  combustion,  in  which  the  combination  of 
preferential  diffusion  and  the  strain-rate  field  leads  to  Tmax  well  above,  or  below, 
Tf.  This  difference  is  clearly  attributable  to  the  hypergohc  nature  of  the  reacting 
system.  Detailed  reactant-flux  analysis  revealed  that  the  ability  to  estabhsh  strain 
rates  well  above  values  encountered  in  the  combustion  of  typical  fuels  can  result 
in  convective  fluxes  that  are  substantially  larger  than  diffusive  fluxes,  everywhere. 
Thus,  convection  is  the  dominant  transport  mechanism  with  all  properties  trans¬ 
ported  at  the  same  rate,  leading  to  a  “locally-conserved”  system. 

The  results  summarized  in  Fig.  20  indicate  that  only  a  negligible  fraction  of  the 
normalized  temperature  rise,  6,  can  be  expected  to  be  contributed  by  combustion  in 
non-premixed  strained  diffusion-flame  regions,  at  strain  rates  of  order  15,000/s,  and 
greater.  Yet,  such,  and  substantially  greater,  values  for  the  strain  rate  are  expected 
in  the  high-speed,  chemically-reacting  shear  layers  discussed  earher.  In  particular, 
the  amount  of  chemical  product  measured  in  these  flows  is  not  reconcilable  with  the 
predicted  chemical  product  formation  in  strained-flame  regions  in  these  simulations, 
considering  the  strain-rate  values  these  flows  can  sustain.  We  are  forced  to  conclude 
that  the  preponderant  fraction  of  chemical-product  formation,  in  these  flows,  occurs 
in  the  premixed  mode,  or,  equivalently,  a  negligible  voliune  fraction  of  mixed  fluid 
can  be  found  in  non-premixed,  strained  diffusion  layers.  Where  such  layers  do  occur, 
the  conditions  for  this  combustion  mode  will  be  created  as  a  result  of  the  substantial 
reactant  leakage,  predicted  to  emanate  from  high  strain-rate  diffusion-flame  regions 
in  the  present  simulations.  Combustion  in  this  mode  can  provide  the  only  candidate 
mechanism  to  explain  the  high  attained  values  of  the  normalized  temperature  rise 
observed  experimentally  (Eq.  14). 

Additional  details  of  this  part  of  the  work  were  documented  in  Egolfopoulos, 
Dimotakis,  and  Bond  (1996). 
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7.  Numerical  simulation  of  flows  with  strong  fronts 

Work  on  the  numerical  simulation  of  flows  with  strong  fronts  continued,  based 
on  the  ideas  and  simulation  technology  in  terms  of  Riemann  Invariant  Manifolds 
(Lappas  1993,  developed  under  a  previous  AFOSR  grant,  Lappas  et  al.  1993,  and 
Lappas  et  al.  1998).  In  particular,  under  sponsorship  of  this  Grant,  this  numerical- 
simulation  methodology  was  extended  to  flow  systems  with  source  terms  in  the 
equations  of  motion,  as  occur  in  chemically-reacting  systems  (Papalexandris  1997), 
for  example.  The  new  developments  permitted  the  reliable  simulation  of  xmsta- 
ble  (underdriven)  one-dimensional  detonation  phenomena,  using  a  Zeldovich-von 
Neuman-Doering,  one-step  chemical-reaction  (ZND)  detonation  model  (Papalexan¬ 
dris  et  al.  1997a),  and,  more  recently,  the  extension  to  two-dimensional  detonations 
(Papalexandris  et  al.  1997b). 

A  key  element  in  these  developments  has  been  the  extension  of  the  classical 
theory  of  characteristics  permitting  the  numerical  simulation  of  various  strong-front 
phenomena  without  resorting  to  operator  dimensional-  as  well  as  flow/chemistry- 
operator  splitting.  By  way  of  example,  this  permitted  the  simulation  of  one-  and 
two-dimensional  detonation  phenomena  with  a  general,  maathematically-consistent 
set  of  equations,  without  resorting  to  added  numerical  viscosity  or  other  numerical- 
dissipation  artifices. 

The  discussion  below  describes  the  application  of  this  methodology  to  unsteady 
detonation  fronts  over  wedge  geometries  in  hypersonic  flow. 


7.1  Numerical  study  of  wedge-induced  detonations 

In  this  section,  numerical  simulations  of  wedge-induced  detonations,  using  the 
proposed  unsplit  scheme,  are  presented.  The  wedge  is  placed  instantly  in  uniform 
flow  of  a  reactive  gas.  A  shock  forms  immediately  at  the  wedge.  If  the  surrounding 
gas  were  inert,  this  shock  would  be  oblique,  at  an  angle  that  is  prescribed  with  re¬ 
spect  to  the  centerline  of  the  wedge.  However,  because  the  gas  is  reactive,  the  shock 
is  curved  due  to  dilatation  of  the  reacting  material  behind  the  shock.  In  particular, 
a  detonation  wave  can  be  established  downstream,  if  the  shock  temperature  is  high 
enough. 

The  wedge  angle,  0,  is  an  important  parameter  in  this  problem.  It  is  expected 
that  for  small  wedge  angles  the  shock  turns  smoothly  and  the  flow  far  downstream 
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consists  of  an  oblique  ZND  wave,  i.e.,  a  ZND  wave  with  a  non-zero  transverse 
velocity  component.  For  small  wedge  angles,  both  the  inert  and  equilibrium  shock- 
polars  admit  solutions  for  the  shock  angle,  The  shock  near  the  tip  is  essentially 
inert  and  its  angle  can  be  computed  from  the  inert  shock-polax.  The  angle  of  the 
ZND  wave,  far  downstream,  can  be  computed  from  the  equilibrium  shock-polar. 
Given  the  state  ahead  of  the  shock,  denoted  by  the  subscript  “1”,  one  can  determine 
the  flow  variables  behind  the  shock,  denoted  by  the  subscript  “s”,  by  employing  the 
standard  kinematic  relations  for  oblique  shocks: 


^  1  +  7^1.  ±  -  1)"  -  2(7  +  l)Mf„  ;o/(cprr) 

(7 + ml 


(15a) 


El  ^  p  = 

Ps  tan  /? 


(15b) 


P.  =  Pi  +  “I  (1  -  F)  , 

where  Mjn  is  the  normal  Mach  number  ahead  of  the  shock: 

sin  /3 


Mi„  = 


n/tT 
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For  wedge  angles  larger  than  a  certain  value  (but  small  enough  so  that  the 
equilibriiun  shock-polar  admits  a  solution  for  l3),  the  shock  carmot  turn  smoothly, 
and  a  strong  explosion  is  expected  to  occur  on  the  front.  This  explosion  is  caused 
by  the  interaction  of  pressure  waves  inside  the  reaction  zone.  These  pressme  waves 
are  emitted  from  the  points  at  which  the  material  near  the  wedge  burns  rapidly  due 
to  thermal  rtmaway. 


7.2  Detonation  induced  by  long  wedges 

Three  different  cases  have  been  examined  numerically.  Schematic  of  the  geo¬ 
metrical  conflguration  and  the  computational  domain  is  given  in  Fig.  24.  The  wedge 
is  assumed  to  be  long  so  that  the  wedge  comer  has  no  effect  on  the  flow-field  near 
the  reaction  zone.  Inflow  conditions  have  been  assigned  at  the  left  boundary  and  at 
the  first  7  cells  of  the  bottom  botmdary.  Reflecting  conditions  have  been  assigned 
at  the  rest  of  the  bottom  boundary,  and  outflow  conditions  have  been  assumed  at 
the  top  and  right  boxmdaries.  The  flow  at  these  boundaries  is  supersonic,  and  one 
must  ensure  that  all  information  for  the  evaluation  of  the  boundary  fluxes  comes 
from  inside  the  domain.  This  is  achieved  by  copying  boimdary-cell  values  to  their 
corresponding  dummy  cells. 


computational  domain--^ 


Fig.  24  Schematic  of  the  computational  domain  for  the  problem  of  a  detonation 
induced  by  a  long  wedge. 

The  values  of  the  specific  heat  ratio,  7,  and  the  heat-release  coefficient,  qo,  are 
set  as  follows: 

7  =  1.2  ,  go  =  50  . 

Upstream,  the  pressure  and  density  of  the  gas  are  set  to  unity.  All  the  variables 
and  parameters  of  the  system  of  equations  have  been  made  dimensionless  with 
respect  to  this  upstream  state.  The  results  of  the  simvJations  are  plotted  with 
the  same  format  that  was  used  in  the  previous  section,  i.e.,  contour  plots  of  the 
pressure,  temperature,  vorticity,  and  reactant  mass  fraction  are  presented.  Each 
plot  consists  of  30  contours,  equally  distributed  between  the  extremal  values,  except 
for  the  plots  of  the  reactant  mass  fraction,  which  contain  11  contours,  at  the  levels 
z  =  0.01 , 0.1 , 0.2, . . . ,  0.9,  0.99.  In  these  simulations,  the  Courant-Priedrichs-Lewy 
(CFL)  number,  U  At /Ax,  where  U  is  the  flow  speed,  and  Ax  and  At  are  the 
magnitudes  of  the  spatial  and  temporal  increments,  is  set  to  0.70. 
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Case  A 

This  is  a  low  activation-energy  case,  *.e., 

Ea  =  10  ,  K  =  3.1245 . 

As  a  first  test,  the  wedge  angle  is  set  to  ^  =  20®.  The  upstream  velocity  of  the  gas 
is  ui  =  12.171.  Given  these  parameters  and  upstream  conditions,  the  theoretical 
prediction  is  that  far  downstream  this  detonation  reduces  to  an  oblique  ZND  wave 
of  overdrive  factor  /  =  1.2  and  at  an  angle  ^  =  34.02°.  The  equivalent  channel-flow 
problem  was  examined  in  the  previous  section  (case  E). 

The  computational  domain  is  resolved  in  960  x  400  cells.  This  resolution  results 
(nominally)  in  8  points  per  half-reaction  length,  i.e.,  the  distance  required  for  half 
the  chemical  reaction  to  be  completed  locally,  for  the  steady  (ZND)  solution.  In  this 
simulation,  the  shock  wave  turns  smoothly  until  it  reaches  the  asymptotic  angle. 
Contour  plots  for  this  case  are  presented  in  Fig.  25a.  These  plots  are  taken  at 
t  =  24.0.  No  change  in  the  flow  variables  is  observed  at  later  times,  which  imphes 
that  the  part  of  the  flow-field  that  is  covered  by  the  computational  domain  reaches 
a  steady  state. 

In  the  area  near  the  tip  of  the  wedge,  no  chemical  reactions  have  had  the  chance 
to  occur  and  the  shock  is  essentially  inert.  The  pressure  and  particle  velocity  in  this 
region  are  almost  constant  along  a  streamline,  and  chemical  reactions  take  place 
due  to  the  phenomenon  of  thermal  runaway.  The  temperature  increase  across  the 
oblique  shock  is  small,  and,  consequently,  the  somce  term  in  the  species  equation 
remains  too  small  to  initiate  a  rapid  reaction.  This  source  term,  however,  is  not 
zero,  so  the  material  does  bum,  if  very  slowly. 

The  region  of  slow  burning  is  the  induction  zone.  At  the  end  of  the  induction 
zone,  the  temperature  has  risen  high  enough  to  initiate  and  sustain  fast  burn¬ 
ing  of  the  gas.  As  a  result  of  the  fast  burning  of  the  material  near  the  wedge, 
pressure  waves  are  transmitted  to  the  shock  front.  These  waves  interact  with  the 
slow-burning  region  behind  the  shock,  causing  fluid  elements  to  bum  faster  and 
dilate,  and  the  shock  front  to  turn.  Recall  that  the  phenomenon  of  thermal  nm- 
away  was  also  encountered  in  the  study  of  one-dimensional  detonations  near  the 
Chapman-Jouguet  (C-J)  point.  This  thermal-runaway  mechanism  is  responsible 
for  the  formation  of  pockets  of  unreacted  material  in  such  flows,  as  discussed  in 
Papalexandris  (1997). 
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Fig.  25a  Case  A.,  9  =  20°:  contoiir  plots  of  the  flow  variables  at  t  =  24.0. 


Fig.  25b  Case  K,0  =  27®:  contoxir  plots  of  the  flow  variables  at  <  =  50.0 
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Next,  the  wedge  angle  is  increased  to  ^  =  27°.  This  is  the  maximum  angle 
for  which  the  equilibrium  shock-polar  admits  a  solution  for  /?.  For  a  detonation 
with  an  overdrive  factor  /  =  1.2  to  occur,  the  upstream  velocity  has  to  be  uj  = 
9.255.  The  asymptotic  limit  of  the  shock  angle  is  readily  found  to  be  =  53.7°. 
The  computational  domain  for  this  simulation  consists  of  1020  x  402  cells,  yielding 
a  nominal  resolution  of  6  points  per  half-reaction  length.  Contour  plots  of  the 
variables  for  this  problem,  at  t  =  50.0,  are  given  in  Fig.  25b.  In  these  plots,  several 
triple  points  can  be  seen  to  have  formed  along  the  leading  front.  These  triple  points 
move  with  different  velocities.  Therefore,  collisions  between  the  triple  points  will 
eventually  occur  (most  likely  outside  the  area  covered  by  the  computational  mesh), 
forming  the  well-known  cellular  structures  behind  detonation  waves.  Formation  of 
triple  points  should  have  been  observed  in  the  previous  simulation  (when  9  =  20°), 
had  the  computational  domain  been  large  enough. 


Case  B 

The  activation  energy  and  the  stiffness  coefficients  are  now  set  at: 

=  50  ,  K  =  99.762  . 

The  wedge  angle  is  0  =  20°.  The  upstream  velocity  is  ui  =  20.58.  Theoretically,  the 
shock  angle  should  tend  asymptotically  to  0  =  27.9°.  In  this  limit,  the  detonation 
is  an  oblique  ZND  wave  with  an  overdrive  factor  of  /  =  2.0. 

Three  different  resolutions  have  been  used  for  this  problem.  In  the  first  test, 
the  computational  domain  consists  of  285  x  60  cells  (Fig.  26a).  In  the  second  test,  it 
consists  of  570  x  120  cells  (Fig.  26b).  In  the  final  test,  it  consists  of  1140  x  240  cells 
(Fig.  26c),  corresponding  to  a  nominal  resolution  of  6  points  per  half-reaction  length 
of  the  one-dimensional,  steady-state  solution.  Again,  the  shock  turns  smoothly  until 
it  reaches  a  steady  state,  very  close  to  the  asymptotic  limit. 

The  results  presented  in  Figs.  26  are  taken  at  t  =  18.0.  By  this  time,  the  solu¬ 
tion  has  already  reached  the  steady  state.  It  is  expected,  however,  that  formation  of 
triple  points  occurs  further  downstream.  Unfortunately,  the  current  constraints  on 
computational  resources  did  not  allow  a  laxge-enough  domain  to  be  implemented  so 
as  to  include  the  region  where  triple-points  occur.  The  same  steady-state  solution 
has  been  computed  with  all  three  meshes.  As  expected,  the  shock  profiles  on  the 
coarse  meshes  are  more  smeared  than  the  ones  on  the  fine  mesh. 
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reoctont  mass  fraction 


Fig.  26a  Case  B,  0  =  20°:  contour  plots  of  the  flow  variables  at  <  =  18.0.  Resolu¬ 
tion,  285  X  60  cells. 

Next,  the  wedge  angle  is  increased  to  ^  =  35°,  which  is  near  the  maximum  angle 
for  which  the  equilibrium  shock-polar  admits  a  solution.  As  before,  three  different 
mesh  sizes  have  been  used  for  this  problem,  consisting  of  560  x  64  cells,  560  x  128 
cells,  and  1140x240  cells,  respectively.  The  finer  computational  domain  corresponds 
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reoctont  moss  fraction 


Fig.  26b  Case  B,  0  =  20®:  contour  plots  of  the  flow  variables  at  <  =  18.0.  Resolu- 
tion,  570  X  120  cells. 


to  a  nominal  resolution  of  8  points  per  half-reaction  length.  The  upstream  velocity 
is  wi  =  11.509.  The  asymptotic  limit  is  a  detonation  with  an  overdrive  factor  of 
/  =  1.2  and  at  an  angle  of  ^  =  56.8°.  Results  for  this  simulation,  at  t  =  50.0, 
are  given  in  Figs.  26d,  26e,  and  26f.  The  flow  cannot  turn  smoothly  in  this  case. 
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Fig.  26c  Case  B,  ^  =  20°:  contour  plots  of  the  flow  variables  at  <  =  18.0.  Resolu¬ 
tion:  1140  X  240  cells. 


because  of  the  high  value  of  the  wedge  angle. 

As  a  result,  a  strong  explosion  takes  place  at  the  front.  The  center  of  the 
explosion  is  a  triple  point.  The  incident  shock  and  the  Mach  stem  are  the  two 


temperature 


FlG.26e  Case  B,  ^  =  35°:  contour  plots  of  the  flow  variables  at  t  =  50.0.  Resolu¬ 
tion,  560  X  128  cells. 

parts  of  the  main  front,  below  and  above  the  triple  point,  respectively.  Another 
(reflected)  shock  emanates  from  the  triple  point,  which  hits  the  wedge  and  reflects 
back.  Additionally,  a  contact  discontinuity  (shear  layer)  is  formed  between  the 
Mach  stem  and  the  reflected  shock.  The  material  behind  the  incident  shock  hums 
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Fig.  26f  Case  B,  ^  =  35°:  contour  plots  of  the  flow  variables  at  t  =  50.0.  Resolution 
1120  X  256  cells. 

due  to  thermal  runaway,  but  the  material  behind  the  Mach  stem  bums  fast  because 
of  the  high  temperature  rise  due  to  the  explosion.  Consequently,  there  are  strong 
density  and  temperature  gradients  across  the  shear  layer.  The  shear  layer  becomes 
unstable  very  quickly  and  generates  strong  vortical  structures  that  are  convected 
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downstream,  as  seen  in  Fig.  26f.  Diffusion  in  these  simulations  is  introduced  by  the 
discretization  of  the  equations  and  the  truncation  error  of  the  scheme.  Consequently, 
the  size  and  the  velocity  of  the  vortical  structures  are  determined  by  the  the  size  of 
the  grid  and  the  implicit  artificial  viscosity  of  the  algorithm.  The  convective  Mach 
number  for  the  shear  layer  lies  between  0.5  and  0.7,  so  the  layer  can  be  considered 
moderately  compressible  by  the  usual  criteria. 

The  flow  between  the  incident  shock  and  the  wedge  is  steady.  Furthermore,  the 
region  between  the  shear  layer  and  the  Mach  stem  is  subsonic  in  the  vicinity  of  the 
triple  point.  A  sonic  line  emanates  from  the  triple  point  and  moves  along  the  shear 
layer.  After  some  distance,  the  sonic  line  turns  upward  and  ends  up  in  the  Mach 
stem.  Beyond  this  pocket,  the  flow  becomes  supersonic  and  should  have  all  the  typ¬ 
ical  characteristics  of  planar  detonations,  such  as  formation  of  colliding  triple  points 
and  transverse  waves,  and  vorticity  generation  at  the  main  front.  Qualitatively,  the 
results  on  the  three  different  meshes  are  the  same.  Interestingly,  the  location  of 
the  explosion  at  the  front  is  the  same  regardless  of  resolution.  The  results  on  the 
coarse  meshes,  however,  are  more  diffuse  than  than  the  ones  on  the  fine  mesh,  as 
expected.  The  shock  profiles  axe  more  smeared,  and  the  shear-layer  structures,  as 
well  as  the  reaction  zone,  are  not  sufficiently  resolved. 

Since  the  equilibrium  shock-polar  admits  a  solution  for  the  shock  angle,  it  is 
expected  that  far  downstream  the  Mach  stem  reduces  to  an  oblique  ZND  wave,  at  an 
angle  /3  =  56.8°.  The  instability  mechanisms  that  were  encountered  in  detonations 
in  channel  flows  (Papalexandris  1997)  and  lead  to  the  formation  of  triple  points 
should  also  appear  on  the  flow-field  of  oblique  detonations.  As  demonstrated  in 
the  previous  section,  two-dimensional  detonations  are  intrinsically  unstable.  This 
region,  however,  is  too  far  downstream  to  be  included  in  the  computational  domain. 

For  higher  wedge  angles,  the  equilibrium  shock  polar  cannot  give  a  solution  for 
13.  In  such  situations,  a  strong  explosion  of  the  main  front  is  also  expected  to  occur. 
Since  there  is  no  solution  to  the  equilibrimn  shock-polar,  the  Mach  stem  cannot 
reduce  to  an  oblique  ZND  wave  downstream,  and  is  everywhere  curved.  For  even 
higher  wedge  angles,  the  shock-polars  do  not  admit  a  solution,  and  the  main  front 
detaches  from  the  wedge.  Such  high  wedge  angles  have  not  been  considered  in  the 
present  work. 
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Case  C 

The  activation  energy  is  the  same  as  in  Case  B,  but  the  stiffness  coefficient  has 
now  been  increased: 

Ea  =  50  ,  K  =  230.75  . 

The  wedge  angle  is  ^  =  20°.  The  upstream  velocity  is  ui  =  18.051.  The  theoretical 
prediction  is  that,  far  downstream,  the  detonation  will  be  a  ZND  wave  with  an 
overdrive  factor  of  /  =  1.60  and  at  an  angle  /3  =  28.5°.  The  equivalent  channel- 
flow  problem  was  examined  in  Case  C  in  Papalexandris  (1997).  The  computational 
domain  consists  of  1200  x  200  cells,  corresponding  to  a  nominal  resolution  of  8  points 
per  half-reaction  length  of  the  one-dimensional,  steady-state  solution.  As  in  case 
B,  the  shock  turns  smoothly  until  it  reaches  a  steady  state.  Contom:  plots  of  the 
flow  variables  are  given  in  Fig.  27a.  These  results  are  taken  at  time  t  =  18.0.  No 
change  in  the  flow  variables  could  be  observed  after  that  time.  The  shock  angle  at 
the  right  boundary  is  very  close  to  the  asymptotic  limit  of  /?  =  28.5°. 

The  case  of  a  higher  wedge  angle,  namely  9  =  30°,  has  also  been  considered. 
The  computational  domain  consists  of  1440  x  240  cells,  corresponding  to  a  nominal 
resolution  of  8  points  per  half-reaction  length.  The  upstream  velocity  is  now  set  at 
ui  =  12.035.  Under  this  initial  condition,  the  theoretical  prediction  is  a  detonation 
of  overdrive  factor  /  =  1.6,  at  a  shock  angle  ^  =  45.7°.  Results  for  this  simulation, 
taken  at  time  t  =  36.0,  are  given  in  Fig.  27b.  These  results  are  similar  to  the  ones 
obtained  in  the  high-angle  simulation  of  Case  B.  The  basic  features  discussed  above, 
z.e.,  the  explosion  on  the  front,  the  development  of  an  unsteady  shear  layer,  and  the 
formation  of  a  subsonic  pocket  behind  the  Mach  stem  in  the  vicinity  of  the  triple 
point,  can  also  be  observed  in  this  simulation. 

Li  et  al.  (1994)  also  presented  simulations  for  wedge-induced  detonations  using 
the  flux-corrected-transport  (FCT)  algorithm,  on  a  domain  of  400  x  150  cells.  Their 
numerical  results  seem  to  agree  qualitatively  with  the  results  obtained  with  the 
proposed  unsplit  algorithm.  In  particular,  they  also  observed  that,  for  small  wedge 
angles,  the  main  front  turns  smoothly,  while  an  explosion  occurs  at  the  front,  if  large 
wedge  angles  are  considered.  More  detailed  comparisons  between  the  two  studies 
cannot  be  made  because  Li  et  al.  considered  a  different  combustion  model  and  did 
not  use  dimensionless  quantities  for  their  simulations. 
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Fig.  27a  Case  C,  ^  =  20°:  contour  plots  of  the  flow  variables  at  t  =  18.0. 

7.3  Detonations  induced  by  short  wedges 

In  the  cases  examined  and  discussed  above,  it  was  assumed  that  the  wedge  is 
long  enough  for  the  explosion  to  take  place  upstream  of  the  comer  of  the  wedge. 
In  such  cases,  the  location  of  the  explosion  and  the  flow-field  in  that  neighborhood 
are  determined  completely  by  the  kinetics  of  the  reaction.  If,  however,  the  wedge  is 
not  long  enough,  then  the  explosion  will  take  place  near  the  comer.  The  effect  of 
the  corner  in  such  situations  has  also  been  studied  numerically,  and  the  results  are 
described  below. 
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vorticity 


reactant  mass  fraction 


Fig.  27b  Case  C,  ^  =  35°:  contour  plots  of  the  flow  variables  at  t  =  36.0. 


The  heat  release  from  the  chemical  reaction  raises  the  temperature  of  the  fluid, 
while  the  expansion  at  the  comer  decreases  it.  These  two  mechanisms  “compete” 
against  each  other.  Furthermore,  including  the  downstream  comer  in  the  compu¬ 
tational  domain  introduces  a  second  characteristic  length,  the  height  of  the  wedge, 
hw,  into  the  dynamics,  in  addition  to  the  half-reaction  length.  The  combination 
of  these  two  length-scales,  with  the  wedge  angle,  6,  determines  which  of  the  two 
mechanisms  will  dominate. 


temoerature 


reactant  mass  fraction 


FlG.26d  Case  B,  ^  =  35®:  contour  plots  of  the  flow  variables  at  i  =  50.0.  Resolu- 
tion,  280  X  64  cells. 

In  the  numerical  simulations,  the  activation  energy  and  the  stiffness  coefficient 
are  set  at  E^.  =  50  and  K  =  99.762,.  respectively.  The  wedge  angle  is  ^  =  35°.  The 
upstream  state  of  the  fluid  is 
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The  flow  produced  by  these  parameters  and  upstream  conditions  in  the  case  of  a 
long  wedge  was  examined  above  (c/.  Fig.  26d). 

Three  wedge  heights  are  considered.  First,  the  wedge  height  is  set  at  hy,  =  70.0. 
The  computational  domain  consists  of  930  x  330  cells.  The  length  of  the  domain 
is  155.0,  and  its  width  from  the  upper  wall  of  the  wedge  is  55.0.  Results  for  this 
simulation,  at  t  =  46.0,  are  shown  in  Fig.  28a.  In  this  case  the  explosion  occurs 
upstream  with  respect  to  the  corner.  The  expansion  at  the  corner  does  not  affect  the 
explosion  because  the  flow  in  the  corner  is  supersonic.  The  leading  shock  is  expected 
to  be  reduced  far  downstream  to  a  C-J  wave,  as  in  one-dimensional  detonations 
initiated  by  a  moving  piston  that  comes  suddenly  to  rest. 


Fig  .  28a  Flow  past  a  wedge  of  height 

of  the  flow  variables  at  t  =  46.0. 


reactant  mass  fraction 


70.0  and  angle  6  =  35°.  Contour  plots 


The  flow-field  in  the  neighborhood  of  the  inert  shock  and  the  explosion  is  the 
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same  as  in  the  corresponding  case  with  a  long  wedge;  see  Fig.  26b.  The  interaction 
of  the  reflected  shock  with  the  shear  layer  causes  the  fluid  to  decelerate,  thus  forming 
another  subsonic  region  in  the  flow-fleld,  besides  the  one  in  the  vicinity  of  the  triple 
point.  The  expansion  that  takes  place  at  the  comer  affects  the  evolution  of  the 
shear  layer.  As  the  fluid  below  the  shear  layer  expands,  the  pressure  and  density 
drop,  generating  even  higher  entropy  gradients  across  the  layer.  The  expansion  at 
the  comer  then  leads  to  an  increased  amount  of  vorticity  generation  across  the  shear 
layer. 

The  wedge  height  was  then  reduced  to  hy,  =  50.0.  For  this  case,  the  compu¬ 
tational  domain  consists  of  810  x  420  cells.  The  length  of  the  domain  is  135.0  and 
its  width,  as  measured  from  the  upper  wall  of  the  wedge,  is  70.0.  Results  from 
this  simulation  are  shown  in  Fig.  28b,  at  time  t  =  46.0.  In  this  case,  the  explosion 
takes  place  downstream  of  the  corner.  The  flow  in  the  subsonic  area  in  the  vicinity 
of  the  triple  point  (between  the  Mach  stem  and  the  shear  layer)  is  influenced  by 
the  expansion  at  the  corner.  Additionally,  the  expansion  at  the  corner  affects  the 
curvature  of  the  leading  front  and  the  reflected  shock. 

A  fluid  element  moving  parallel  to  the  wedge  does  not  have  time  to  increase 
its  temperature  substantially,  via  the  thermal-runaway  mechanism,  because  of  the 
small  length  of  the  wedge.  Consequently,  it  remains  almost  unreacted  when  it 
reaches  the  head  of  the  expansion.  The  expansion  produces  a  further  decrease  of 
the  temperature,  which  delays  the  initiation  of  the  reaction  even  more.  Material  on 
the  upper  wall  of  the  wedge  has  remained  only  partially  reacted  (about  20%),  as  a 
result  of  the  temperature  decrease  caused  by  the  expansion  at  the  corner. 

As  a  final  test,  the  wedge  height  is  decreased  further  to  hw  =  25.0.  The 
simulation  is  performed  on  a  domain  of  1266  x  270  cells.  The  length  of  the  domain 
is  210.0  and  its  width  is  45.0.  Results  from  this  simulation,  taken  at  t  =  10.0, 
are  presented  in  Fig.  28c.  The  expansion  at  the  corner  reduces  the  temperature 
so  much  that  the  gas  near  the  wedge  remains  almost  imreacted  because  the  time 
needed  for  rapid  reaction  via  thermal  runaway  becomes  very  large.  As  a  result, 
a  detonation  cannot  be  established  and  the  shock  wave  is  expected  to  be  reduced 
to  a  Mach  wave,  downstream.  There  is  also  a  small  pocket  of  slightly-reacted 
material  near  the  wedge.  Results  from  early  times  indicate  that  the  formation 
of  this  pocket  is  a  transient  phenomenon  caused  by  the  interaction  of  the  shock 
wave  and  the  expansion  at  the  comer,  in  the  beginning  of  the  simulation.  The 
pocket  is  convected  downstream  with  the  fluid  velocity.  The  TninimiiTn  value  of 
the  reactant  mass  fraction  inside  the  pocket  is  ^  ~  0.9.  If  the  material  reacted 
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Fig.  28b  Flow  past  a  wedge  of  height  h^,  =  50.0  and  angle  6  =  35°.  Contour  plots 
of  the  flow  variables  at  t  —  46.0. 


completely,  pressure  waves  would  be  transmitted  to  the  shock  wave,  and  a  C-J 
detonation  could  be  established.  The  reaction  rate,  i  =  —Kzexp{—Ea,/T)^  inside 
the  pocket  decays  with  time,  however,  indicating  that  the  reaction  process  will  not 
be  completed  there.  After  some  time,  the  pocket  exits  the  computational  domain 
and  the  shock  front  assumes  a  fixed  position.  No  change  in  the  flow  variables  are 
subsequently  observed.  The  material  all  along  the  wedge  will  remain  only  partially 
reacted  and  a  detonation  will  not  be  estabhshed. 
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Fig.  28c  Flow  past  a  wedge  of  height  =  25.0  and  angle  ^  =  35°.  Contour  plots 
of  the  flow  variables  at  t  =  10.0. 


7.4  Concluding  remarks 

The  proposed  unsplit  numerical  scheme  has  been  used  for  the  numerical  study 
of  two-dimensional  detonating  flows.  The  scheme  is  based  on  a  decomposition  of  the 
governing  equations  that  is  mathematically  consistent  and  appears  capable  of  cap¬ 
turing  important  details  of  the  structure  of  the  resulting  flow-fields,  some  of  which 
have  not  been  easy  to  simulate  in  the  past.  No  explicit  artificial  viscosity  mech¬ 
anisms,  or  any  other  of  the  usual  “Axes”  (such  as  entropy-fixes  or  flux-splitting), 
have  been  employed. 

With  this  scheme,  it  was  verified  that  two-dimensional  detonations  propagat¬ 
ing  in  narrow  channels  are  intrinsically  unstable  and  exhibit  chaotic  behavior  (Pa- 
palexandns  1997).  They  are  characterized  by  the  presence  of  cellular  patterns. 
These  patterns  are  formed  by  the  transverse  waves  of  the  triple  points  of  the  main 
front.  The  slip  lines  emanating  from  the  triple  points  are  rolled-up  vortex  sheets. 
They  detach  from  the  leading  shock  when  triple-point  collisions  occtur.  Vortical 
structures  are  also  generated  by  shock  curvature,  but  their  strength  decreases  sub¬ 
stantially  within  the  chemical  reaction  zone. 

In  addition  to  a  host  of  detonation  phenomena  (Papalexandris  1997;  Papalexan- 
dris,  Dimotakis,  and  Leonard  1997a,b)  a  numerical  study  of  shock-wedge  induced 
detonations  was  also  performed,  as  described  above.  These  simulations  demon¬ 
strated  that,  for  small  wedge  angles,  the  shock  attached  to  the  wedge  truns  smoothly 
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to  an  oblique  ZND  wave.  For  high  wedge  angles,  however,  such  a  smooth  tmn  is 
not  possible  and  an  explosion  takes  place  at  the  front.  The  center  of  the  explosion 
is  a  stationary  triple  point.  A  shear  layer  emanates  from  the  triple  point,  that  is 
hydrodynamicaJly  unstable. 

The  effect  of  the  downstream  corner  of  the  wedge  was  also  studied  ntimerically. 
It  was  shown  that  when  the  explosion  of  the  leading  front  takes  place  upstream  of 
the  comer,  the  expansion  at  the  corner  does  not  affect  the  evolution  of  the  front, 
which  reduces  to  a  C-J  wave.  When  the  explosion  occturs  downstream  of  the  corner, 
the  curvature  of  the  front  and  the  reaction  process  depend  on  the  expansion  at  the 
comer.  It  appears  that  for  wedge  heights  small  enough,  a  detonation  cannot  be 
established  downstream,  and  the  front  decays  downstream  to  a  Mach  wave. 
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